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LARGE EDDY SIMULATION OF INCOMPRESSIBLE TURBULENT CHANNEL FLOW 


Abstract 

The three-dimensional, time-dependent primitive equations of motion 
have been numerically integrated for the case of turbulent channel flow. 

For this purpose, a partially implicit numerical method has been devel- 
oped. An important feature of this scheme is that the equation of con- 
tinuity is solved directly. The residual field motions were simulated 
through an eddy viscosity model, whereas the large-scale field was ob- 
tained directly from the solution of the governing equations. 16 uniform 
grid points were used in each of the streamwise and spanwise directions, 
and 65 grid points with non-uniform spacings in the direction normal to 
the walls. An important portion of the initial velocity field was ob- 
tained from the solution of the linearized Navier-Stokes equations. The 
pseudospectral method was used for numerical differentiation in the hori- 
zontal directions, and second-order finite-difference schemes were used 
in the direction normal to the walls. 

It has been shown that the Large Eddy Simulation technique is capable 
of reproducing some of the important features of wall-bounded turbulent 
flows. The overall agreement of the computed mean velocity profile and 
turbulence statistics with experimental data is satisfactory. The resolv- 
able portions of the root-mean square wall pressure fluctuations, pressure 
velocity-gradient correlations, and velocity pressure-gradient correlations 
are documented. 
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Chapter I 
INTRODUCTION 


1.1 Historical Background 

It has been known for some time that any turbulent flow contains 
structures ("eddies") in a wide range of spatial as well as temporal 
scales. It is also generally recognized that large eddies differ 
markedly from one flow type to another (e.g., jets vs. boundary layers), 
while the small eddies are quite similar in all flows. 

Unfortunately, in the numerical simulation of (high Reynolds number) 
turbulent flows, we find that due to computer limitations one cannot 
resolve all the scales. It is this deficiency which provides the primary 
inducement for the utilization of the large eddy simulation (L.R.S.) 
approach. 

The foundation on which this approach relies concerns the contrast 
between large and small eddy modeling. More specifically, one finds that 
large eddies cannot and should not be modeled, whereas with small eddies 
successful modeling is possible. 

The large eddy simulation method is initiated by the introduction of 
a procedure which separates the small and large scale structures. The 
large scale structures will then be computed explicitly, while the small 
scales are necessarily modeled. 

The problem of decay of homogeneous isotropic turbulence has been 
the subject of extensive study at Stanford University (Kwak et al. 

(1975), Shaanan et al. (1975), Mansour et al. (1977), Ferziger et al. 
(1977)). These studies have shown that with the use of algebraic models 
and a relatively small number of mesh points (16 x 16 x 16 or 
32 x 32 x 32) , homogeneous turbulent flows can be simulated reasonably 
well. 

The first application of the method to problems of engineering 
interest was made by Deardorff (1970) who treated the channel flow 
problem. In his pioneering work, Deardorff showed that a three dimen- 
sional numerical simulation of turbulence is feasible. He was able to 


1 



predict some of the features of turbulent channel flow with a fair amount 
of success. However, as will be clear in the next section, neither 
Deardorff nor the followup work of Schumman (1973) treated the most 
important part of_the flow, namel-y the region very near the wall. It is 
in this region that virtually all of the turbulent energy production 
occurs. By introducing artificial boundary conditions, they, in effect, 
modeled the turbulence production mechanism in the wall region. 

Finally, we note that, concurrent with the present work, Mansour 
et al. (1978) simulated a time developing turbulent mixing layer. They 
showed that essentially all the features of a turbulent mixing layer can 
be reproduced using the L.E.S. approach. 

1.2 Experimental Background 

Many early studies of the structure -of turbulence consisted of 
measurments of the root-mean square and spectra of the turbulent velocity 
fluctuations. Among the measurements that were primarily concerned with 
turbulent boundary layers were those of Townsend (1951) , Klebanoff 
(1954) , Willmarth and Wooldridge (1963) , and for flow near the wall (in 
a pipe) Lauf er (1954) . 

Willmarth made a single, unpublished attempt, in 1960, to bring 
together the then existing results of turbulence-intensity profiles 
of the bound ary l ayer-.o n a s ingle plot (see Willmarth, 1975). The 
curves of Vu ^/u , V v ' ^/u^ , , and V w r 2 / u^, as a function of 
y^/6 (or y + = y w u^_/v) did not agree very well (not within 50%). 

Here, y is the distance to the .wall, u is the shear velocity, 
and 6 is the boundary layer thickness. .Part of the lack of agree- 
ment was attributed to freestream disturbances or differences in the 
methods used to trip the boundary layers. However, in spite of the 
differences between various measurements of turbulence intensity, it 
i s de finitel y est ablish ed th at within a turbulent boundary layer, 

V u ,2 /U > > w' 2 /U > Vv ,2 /U . These differences between the root- 

CO oo 00 

mean-square velocity fluctuations becom e lar ger as one__a pproaches 
the wall. Furthermore, the profiles ^u' 2 and "Vw' 2 have pro- 
nounced local maxima very near the wall. 
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From the measured distributions of turbulence kinetic energy, 
turbulence shear stress, and dissipation, it is possible to obtain a 
turbulence energy balance. Townsend (1951) and Laufer (1954) (among 
others) made such a balance in a boundary layer and pipe flow respec- 
tively. From these data, it can be seen that the production and dissipa- 
tion terms are nearly equal but opposite to each other, and so are the 
terms representing diffusion by turbulence of kinetic energy and of 
pressure energy. Furthermore, it may be noted that the turbulence 
kinetic energy, its production and its dissipation, ail show sharp maxima 
in the buffer region (y+ - 10) near the wall. On the basis of energy 
measurements, Townsend (1956) proposed a two-layer model for the energy 
transformation process. According to this model, the whole layer is 
arbitrarily divided into two parts: (i) an inner layer which is nearly 

in energy equilibrium but within which most of the turbulence production 
takes place, and (ii) an outer layer whose Reynolds stresses retard the 
mean flow but whose principal source of turbulent energy comes from the 
inner layer. 

The level of turbulent intensity in the outer two-thirds of the 
flow is maintained by transport of energy from the inner region since the 
production of energy in the outer region is too small to. balance the 
viscous dissipation and transport losses. Townsend concluded that the 
interaction between the inner and outer layers of the flow may be con- 
sidered as two distinct processes: (i) the transfer of mean-flow energy 

from the outer region to the inner layer at a rate controlled by the 
gradient of Reynolds stresses in the outer layer, and (ii) the transport 
of turbulent energy from the inner layer to the outer layer . 

To gain insight into the mechanics of turbulence production a 
thorough study of the structure of the inner layer was required. 
Runstadler et al. (1959) , (1963) advanced a model for the inner layer 
based, on visual observations using dye and hydrogen bubbles. Their 
studies revealed new features of turbulent boundary layers. In partic- 
ular, they demonstrated that the wall layer is not two dimensional and 
steady; rather it consists of relatively coherent structures of low and 
high speed streaks alternating in the spanwise direction over the entire 
wall. The non-dimensional mean spacing between the low speed streaks 
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was shown to have a universal correlation for fully turbulent layers 
based on wall layer parameters; this is given by the relation 
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The streak pattern is not stationary in space. It migrates and displays 
strong intermittent motion. These intermittent motions involve primar- 
ily the movement of low speed streaks away from the wall. When the 
streak has reached a point corresponding to y+ <_ 8-12 , it begins to 
oscillate. The oscillation grows in amplitude and it is followed by 
breakup. The region where most of the low speed streak breakups are 
observed to occur, i.e., the inner edge of the buffer zone, is the 
region where a sharp peak is seen to occur in the production curve 
(Klebanoff 1954). Kline et al. (1967) and Clark and Markland (1970) 
observed U shaped vortices occasionally in the inner region. In the 
studies of Clark and Markland, an average spanwise spacing of these U 
shaped vortices of A^ - 100 and streamwise spacing of A^ of 440 was 
found. 

Kim et al. (1971) studied bursts using motion pictures of the tra- 
jectories of hydrogen bubbles. From their analysis, they concluded that 
in the region 0 < y + < 100 essentially all the turbulence production 
occurs during bursting. They also observed that during gradual liftup 
of low-speed streaks from the sub-layer, unstable ( inflectional) instan- 
taneous velocity profiles were formed. One of the important findings of 
Kim et al. was that, while the bursting process indeed contributes to 
the turbulent energy, its main effect is to provide turbulence with u’ 
and v* in proper phase to give large positive Reynolds stresses as 
required for the increase in production. 

The findings of Kline and his colleagues were largely confirmed and 
supplemented by the visual studies of Corino and Brodkey (1969) . One of 
their observations was that, after formation of low speed streak a much 
larger high speed bulk of fluid came into view and by "interaction" 
began to accelerate the low speed fluid. The entering high speed fluid 
carried away the slow moving fluid remaining from the ejection process; 
this they called the "sweep" event. 
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The above experimental investigations of the structure of turbulent 
boundary layers are by no means the only ones reported. The number of 
publications on the subject is already very large. Among these is the 
work of Narahari, Rao, Narasimha, and Badri Narayanan (1971), where the 
frequency of occurrence of bursts was studied. Their investigation 
showed that the mean bursting frequency scaled with the outer rather 
than inner flow variables. This was also reported by Kim et al. (1971). 
The recent experimental investigation of Blackwelder and Kaplan (1976) 
studied the near wall structure of the turbulent boundary layer using 
hot-wire rakes and conditional sampling techniques. Among their find- 
ings was that, the normal velocity is directed outwards in the regions 
of strong streamwise-momentum deficit (with respect to the mean velocity) , 
and inwards in the regions of streamwise-momentum excess. This was also 
reported by Grass (1971) . For further details and description of other 
works on the structure of turbulent boundary layers the reader is 
referred to the review articles of Willmarth (1975) and Laufer (1975). 

An entire meeting was recently devoted to review of the state of knowl- 
edge in this area (Abbott 1978) . 

1.3 Motivation and Objectives 

The present study is one in a systematic program investigating large 
eddy simulation of turbulence. In order to extend the available tec.h- 
noj ogy of the L.E.S. approach to wall-bound flows, we chose to study 
incompressible turbulent channel flow. Due to the simplicity of its 
geometry and some experimental advantages, channel flow has been a par- 
ticularly attractive reference flow for both theoretical and experimental 
investigations. As a result, there is a considerable amount of experi- 
mental as well as theoretical findings available for a detailed evalua- 
tion of the large eddy simulation technique. In addition, this flow 
possesses important features of the flows of practical interest. This, 
in turn, allows the evaluation of the L.E.S. approach from a practical 
point of view. 

The specific objectives of this work may be stated as follows: 


5 



a) To develop a numerical method for long time integration of the 
three-dimensional governing equations for the large scale field 
in a turbulent channel flow; 

b) To carry out numerical solution of these equations using a 
simple subgrid scale model; 

c) To evaluate the performance of the Large Eddy Simulation tech- 
nique in reproducing some of the laboratory observations and 
measurements described above, and to compute quantities such as 
pressure velocity gradient correlations that cannot be measured. 

!■* 4 Summary 

The contributions of the present work include: 

a) Demonstration of the inherent numerical problems associated with 
the explicit numerical solution of the dynamical equations of 
motion in primitive form. 

b) Derivation of consistency conditions for the initial velocity 
field such that the Neumann and Dirichlet problems for the pres- 
sure have the same solution. 

c) Development of a new semi-implicit numerical scheme for the 
solution of dynamical equations in primitive form. 

d) Development and use of a new subgrid model in the wall region 
of the turbulent flow. 

e) Development and use of a solution of the Orr-Sommerfeld equation 
for a three-dimensional disturbance as an important part of the 
initial velocity field. 

f) Demonstration that the Large Eddy Simulation technique is 
capable of reproducing many of the important features of the 
turbulent boundary layer. 
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Chapter II 


MATHEMATICAL FOUNDATIONS 

2*1 Definition of Filtered and Residual Fields 

In the large eddy simulation approach, the first and most 
fundamental step is defining the large-scale field. To accomplish this 
task, each author has adopted a slightly different approach, but they 
can be treated within a single conceptual framework as shown by Leonard 
(1974). If f is some flow variable, we decompose it as follows: 

f = f + f * (2.1) 

where f is the large-scale component and f' is the residual field. 
Leonard defined the large scale field as: 

f(x) « f l(x-x’) f(x’) dx' (2.2) 

where G(x-x') is a filter function with a characteristic length A , 
and the integral is extended over the whole flow field. It is to be 
noted that the above form of G (a function of (x-x 1 ) ) is most 
suited for filtering in the directions in which the flow is homogeneous. 
In other words, we point out that the filter function need be neither 
isotropic nor homogeneous and there are many flows (or directions in a 
given flow) in which neither of these properties are desirable. In the 
present work we use the Gaussian filter, 

n ^ 

G(x-x’) = Ufa) exp jJ-6(x i -x|) 2 /A^j (2.3) 

where A_^ = 21k , Ik is the mesh size in the i-direction, and n = 1, 
2, or 3, is the number of dimensions in which the flow is homogeneous. 
Thus in the simulation of the decay of homogeneous isotropic turbulence, 
n = 3 , while in the simulation of turbulent channel flow, we have used 
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n = 2 . A convenient property of a homogeneous filter, G(x-x/) , is 
its commutivxty with partial differentiation operators; using integra- 
tion by parts one can show (Kwak et al. (1975)) : 


df _ 3f 

3x. 9x. 

1 i 


(2.4) 


Due to variation of the physical length scale of turbulence in the 
direction in which the flow is homogeneous, one should not use homoge- 
neous filters in that direction. This is particularly true in turbulent 
boundary layers. Instead, one should use a filter with variable width 
A(r) , where jc is the direction in which the flow is inhomogeneous. 

On the other hand, using a filter with variable width causes some mathe- 
matical difficulties; in particular (2.4) will no longer hold. In 
Appendix A, we explore filters with nonuniform width in some detail. 

Finally, we note that, in the numerical simulation of turbulent 
channel flow, we filter only in the directions in which the flow is 
homogeneous, (streamwise and spanwise directions) i.e., we do not formal- 
ly filter in the direction perpendicular to the walls. The justifica- 
tion for this choice is twofold: 

a) We are using a second order finite difference scheme to 
approximate partial derivatives in the inhomogeneous direction, 
and finite difference shcemes in general have inherent filter- 
ing effect. 

b) The Leonard term is fairly well represented by the truncation 
error of the second order central differencing scheme. (See 
Shaanan (1975)). 

The main disadvantage of this choice is that we do not have a formal 
closed mathematical expression relating the filtered to the unfil- 
tered field. 


2. 2 Dynamical Equations in Primitive Form 

Now let us derive the primitive dynamical equations for the large- 
scale flow field. Starting with the incompressible Navier-Stokes equations. 
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3u. ~ 

1 , 3 

___ + _ — u .u . = 

3t 3x. i 3 


-1 rx 3^u. 

- 1J2- + v i_ 

p 3x. 3x.3x. 

i 3 3 


we can apply the operation ( 2 . 2 ) to get the dynamical equations of large 
scale field. 


9u i „ 3 =“=- 1 3? 3 

-KT" + X U.U. x— - -r T. . + V 

ot 3x. l ] p 3x. 3x . xj 

3 i J 

where we have decomposed as in ( 2 . 1 ) and: 



3x. 3x,. 
3 3 


(2.5) 


T. . = R. . - R, , 5. ./ 3 

xj xj xk xj 

g = P/P + \ k /3 

R. . = u!u! + u!u. + u!u. 

3-J x 3 3 1 i j 

The T„ represents the (negative) subgrid scale stresses and must be 
modeled. We can write (2.5) in the following equivalent form: 


3u. /3u, 3u.\ 


where 


_3P. _ _3_ , U i 

3x. 3x. T ij V 3x.3x. 

i 3 3 3 


( 2 . 6 ) 


P 


+ § 


The rationale for using this form of the equation will be explained 
in Section 3.5. 

In order to calculate the second term on the left-hand side of (2.6), 
we use ( 2 . 2 ) to write: 
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_ /3u. 3u.\ /*+°° 

u j (se - ■ J G< ^ ,) u j 

■3 ^ _co 




dx* 


Note that, here, the filtering and the corresponding integration is 
performed only in the directions in which the flow is homogeneous. Let 
us Fourier transform the above equation (in the homogeneous directions) 
to get: 


A 




A 



(2.7) 


where denotes a Fourier-transformed quantity; a A over a bracket 
means the transform of the bracketed quantity. Thus, given a velocity 
field, u^ , one can compute the term in the brackets on the right-hand 

A 

side of the above equation, Fourier- transform it, multiply it by G , 
and invert the transform to obtain the desired term. 


2.3 Residual Stress Model 

An eddy viscosity model is used for t. . : 

13 


T. . 

13 



S. . 
13 


( 2 . 8 ) 


where 



is the strain rate tensor and is an eddy viscosity associated with 

the residual field motions. In the remainder of this section, we 
present the models used for . Throughout, we assume that the sub- 
grid scale production and dissipation of turbulent kinetic energy are 
equal . 
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Production of the subgrid scale turbulent kinetic energy is given 


by: 



2v m S . . 
T ij 


S. . 
11 


(2.9) 


Inclusion of the experimental observation that, remote from the wall, 

dissipation is controlled only by the largest subgrid-scale eddy param- 

2 

eters such that D = D(q ,£) , coupled with dimensional analysis, 

3 

produce the result first found by Kolmogorov in 1942 that Da q /& . 
Here, q and £ are the characteristic velocity and length scale of 
subgrid scale eddies respectively. Using Prandtl's assumption for eddy 
viscosity, = C^q& , and equating the subgrid production and dissipa- 

tion, we get: 


— — 2 

2C,q£ S.. S.. = q /£ (2.10) 

1 id ij 

From (2.10), we readily obtain: 

q - C 0 £ J2S, J. . 

3 il n 

Again, using Prandtl's assumption, we get: 

V., = (C £) 2 /2S . .S . . (2.11) 

T s ij ij 

This is Smagorinsky ’ s (1963) model, and is to be used in the regions 
away from the solid boundaries. 

On the other hand, very near the wall, the size of the eddies is 

inhibited, and the eddies are of such a size that viscosity can be a 

» 

dissipative agent for the largest eddies. In fact, at the wall, the 

eddy viscosity as well as its gradients should vanish. Under such con- 

2 

ditions viscosity is a factor and D = D(v,q ,£) . Application of 

dimensional analysis to this condition produces the result that 
2 2 

D “ (Vq /£ )f(q£/v). Moreover, at the wall the subgrid scale dissipa- 
tion is given by: 
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D = V 



- - - - - 2 '2 

Thus, in the vicinity of the wall, we assume that D^Vq /£ . Equating 
subgrid scale production and dissipation, we obtain for the inner region 
of the boundary layer: 


V T = (C 2 £ 4 /v) (2S S ) (2.12) 

where C is a constant. 

In order to determine the value of C n , we assume that C , 

2 s 

Smagorinsky 1 s constant, is known from some other calculation e.g., 

simulation of the decay of isotropic turbulence. Strictly speaking, 

there is no rigorous justification that the constant obtained from the 

simulation of a totally homogeneous flow is applicable in the simulation 

» 

of a wall-bounded turbulence with mean shear. Furthermore, in order to 
determine the value of C 2 , several known characteristics of turbulent 
boundary layers will be applied. Among these characteristics is that, 
in the logarithmic section of the layer, the slope of the mean velocity 
profile in the semilogarithmic wall coordinates is 1 /k , where ic is 
the von Karman constant. Hence, in what follows, we give only a rough 
estimate of the value of C 2 , which will be used throughout our simula- 
tion of turbulent channel flow. 

At the edge of the logarithmic section of the boundary layer, (say 
y + =27) , we demand that the inner and outer layer models have the same 
planar mean value. If we nondimens ionalize all the velocities by the 
shear velocity, u , and the lengths by the channel half width, S , 
we have in the logarithmic region: 


/ 


2S. .S. . 
3-J 


9U „ 1 

3y Ky 
w 


(2.13) 


where y is the distance to the lower wall (the lower wall is located 
at y = -1 and the upper wall at y = +1 ) . Note that here, we have 
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assumed that the mean velocity gradient is much larger than all the other 

+ 

velocity gradients- Equating the two models at y = 27 , we obtain: 


C 


2 



(2.14) 


where we have assumed that £ = Ky . Thus , the actual model used for 

J w 

the eddy viscosity at each time step in the calculation is: 


C 2 Re T « 4 (2S ljSlj ) y < ^ 


V t" 


(C £) 2 As. .S. . y > y 

s 13 13 J c 


(2.15) 


Here y c is the coordinate of the first computational grid point away 
from the wall at which the planar average of the two models are closest 
to each other. It is to be noted that, y^ can vary in time and in 
general it does. The same relation as (2.15) is used in the upper half 
of channel (0 <_ y _< 1) . Finally, we turn our attention to the specifi- 
cation of £ . 

Due to the no-slip boundary condition, £ must vanish at the walls. 
Furthermore, due to lack of spatial resolution in the homogeneous direc- 
tions (see Section 3.1), and with no further reasoning, we have used the 
following expression for £ in the simulation of turbulent channel flow: 


£ 



where £' is the Prandtl’s mixing length: 


(2.16) 


0.1 

-< 

V 

.1/k 


w 


^w 

y 5 

w 

1 — i 

« 
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and are the nondimensionalized filter widths in stream- 

wise and spanwise directions respectively, and h 2 is the local grid 
size in the vertical direction. Two remarks are in order. First, due 
to -the particular grid sizes chosen (see Section 3.1), we have the 
following global inequalities: 


h 2 (y) < .1 


Ai > A3 > .1 


(Note that all the lengths are nondimensionalized with respect to 

channel half width 6 ) . Second, we should mention that the expression 

(2.16) for Z is strictly speaking, based on ad hoc foundations and 

more work in this area is strongly recommended (see Chapter V). This 

expression was chosen initially on a trial basis; nevertheless, we did 

not find any alteration of it necessary. Thus, we emphasize that in 

obtaining the computational results presented here, no fine adjustments 

’> of either C g or Z were made. In spite of this, the numerical 

results (see Chapter IV) are satisfactory. It is believed however, that 

an optimum choice for C and Z would somewhat improve the quantita- 

s 

tive results. 


2 .4 Governing Equations for the Large Scale Field 

In the numerical simulation of turbulent channel flow, all the 
variables are nondimensionalized by turbulent shear velocity, u T , and 
the channel half width, 5 . In this case, we solve the following 
equations numerically: 





3P 

3x. 

1 


+ 5 

il 



(2v S .) 
^ T ir 




3x.3x. 
3 3 


(2.19) 


14 



and 


3u . 

3 ^ - 0 <2 ' 20 > 

1 

where Re^ is the Reynolds number based on shear velocity, , and 

channel half width, 6 . Note that the second term on the right-hand 
side of equation (2.19) is the mean pressure gradient imposed on the 
flow. 
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Chapter III 


NUMERICAL METHODS 


3.1 Grid Selection 

For a given number of grid points, N , one has to choose the grid 
size(s) based on the physical properties of the problem at hand. In the 
simulation of the decay of homogeneous isotropic turbulence, for 
example, it is desirable to select the grid size, h , such that the 
filtered field contains as much of the turbulence energy as possible 
(Kwak et al. , 1975) . On the other hand, the length of the side(s) of 
the computational box in the direction (s) in which periodic boundary 
conditions are used should be long enough to include the important large 
eddies (Ferziger et al., 1977). 

In the grid size selection process for the numerical simulation of 
turbulent channel flow, one has to consider the average spanwise and 
streamwise spacing of the turbulent structures in the vicinity of the 
wall (see Section 1.3) as well as the integral scales of turbulence. In 
addition, quantities such as the thickness of the viscous sublayer should 
be taken into consideration. With this in mind we proceed to specify 
our grid system: 

In the vertical direction (-1 <_ y <_ 1) , a nonuniform grid spacing 
is used. The following transformation gives the location of grid points 
in the vertical direction (Mehta, 1977). 


y. = — tanh 

3 a 


E . tanh "^(a) 

L 3 


(3.1) 


where 


Z. = -1 + 2(j-2)/(N-3) j=l,2,...,N 

and N is the total number of grid points in the y direction. Here, a 
is the adjustable parameter of the transformation (0 < a < 1) ; a 
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large value of a distributes more points near the boundary. In our 
computation we have used a = .98346 , and N = 65 . Table 3.1 shows 
the distribution of the grid points in the vertical direction with the 
corresponding values of y+ = y^u /v . Note that in reference to the 
vertical direction, index (or subscript) 1 and N refer to grid points 
just outside the lower and upper walls respectively. 

For the grid selection in the streamwise, x , and spanwise, z , 
directions, one needs to consider the experimentally measured two point 
correlation functions 

R ii ( r ,°, 0 ) = < u^(x,y,z) u ^(x+r,y,z) > 

and 

R ii^°’°’ r ^ = < u i< X ’ y ’ z) u i (x,y,z+r) > 

Here < > denotes the average over an ensemble of experiments. 

The use of periodic boundary conditions in a given direction can be 

justified if the length of the side of the computational box in that 

direction is at least twice the distance r , at which the appropriate 

R. . vanishes, 
n 

Experimental data of Comte-Bellot (1963) , indicates that 

X ± - 6.46 

and 

X 3 = 3.26 

where and X^ are twice the distance, r , beyond which 

R 1;L (r,0,0) and R 1;L (0,0,r) respectively, are negligible. Here 6 
is the channel half width. 

For a complete simulation of the important large scale field , one 
has to select the number of grid points in the streamwise, x , and 
spanwise z , directions with careful consideration to laboratory 
observations. We assume that L and L , the lengths of the computa- 
tional box in the streamwise and spanwise directions, are fixed in 
accordance with the above considerations. As was mentioned in Chapter I, 
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Table 3.1 


GRID DISTRIBUTION, IN THE VERTICAL, y , DIRECTION 


y+ * 


y w = |i+y| 


1 

-1.002 

.002 


2 

-1.000 

.000 

0.000 

3 

- .997219 

.00278 

1.78 

4 

- .993983 

.00602 

3.85 

5 

- .99022 

.00978 

6.26 

6 

- .985847 

.01415 

9.06 

7 

- .980767 

.01923 

12.31 

8 

- .974871 

.02513 

16.09 

9 

- .968035 

.03197 

20.47 

10 

- .960117 

.03988 

25.53 

11 

- .950956 

.04904 

31.40 

12 

- .940372 

.05963 

38.18 

13 

- .928164 

.07184 

45.99 

14 

- .914109 

.08589 

54.99 

15 

- .898 

,102 

65.33 

16 

- .879 

.121 

77.47 

17 

- .858 

.142 

90.91 

18 

- .834 

.166 

106.28 

19 

f 

• 

CO 

o 

.193 

123.57 

20 

- .776 

.224 

143.42 

21 

- .741 

.259 

165.82 

22 

- .702 

.298 

190.79 

23 

- .659 

.341 

218.32 

24 

- .611 

.389 

249.06 

25 

- .559 

.441 

282.35 

26 

- .502 

.498 

318.84 

27 

- .440 

.560 

358.54 

28 

- .374 

.626 

400.80 

29 

- .304 

.696 

445.61 

30 

- .231 

.769 

492.35 

31 

- .156 

.844 

540.37 

32 

- .078 

.922 

590.31 

33 

.0 

1.000 

640.25 


Re^ = 640.25. 
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experimental data indicate that the average (spanwise) streak spacing 
corresponds approximately to A^ =; 100 and the average streamwise 
spacing of the U shaped vortices corresponds to A^ - 440 . Therefore, 
for the channel flow under consideration (see Chapter IV) , the average 
dimensionless distance between the spanwise and streamwise structures 
are: 


^ = 100/Re = 0.156 

0 T 

and 

Ax 

= 440/Re = 0.687 

O T 

respectively. Here Re^ is the Reynolds number based on shear velocity, 
u and channel half width, <5 and is 640 in our simulation. 

Using the above values of and , and assuming that, at 

least four grid points are needed to resolve one wavelength (structure), 
we arrive at the following requirements for the number of grid points in 
x and z directions: 


N = 37 

x 

N = 82 

z 

It is emphasized that the above values for N and N are based on 

X z 

ensemble averaged spacing of the structures. Hence for an adequate 
simulation of the important large scales, the following values for N 
and N z are recommended (with due consideration to the capability of 
present computers): 

N = 32 

x 

N = 128 

z 

In the present numerical simulation of turbulent channel flow, we 
have chosen the following values for the nondimensionalized streamwise 
and spanwise computational box lengths: 
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L = 2ir 
x 


L 

z 


4 

3 


71 


The value of L = -r 7T is somewhat bigger than the above value for 

Z j 

X^/6 . This choice was made with due consideration to stability and 

resolution requirements of linear hydrodynamic stability theory (see 

Section 4.3). In addition, due to computer cost and storage limitations, 

we have used 16 grid points with uniform spacing, in each of the stream- 

wise and spanwise directions. Therefore, the actual grid spacing used 

+ + 

in these directions corresponds to h^ = 251 and = 168 respec- 
tively. Hence, it is clear that we have inadequate resolution, partic- 
ularly in the spanwise direction. 


3.2 Numerical Differentiation 


In the vertical direction, central differencing is employed with 


y j ' y j-l 


and 


variable grid spacing y^. +1 = y^ + h^ +1 where tu 
j = 1,2,...,N (see Section 3.1). The partial derivatives for this case 
are the following expressions with the first truncation error term 
included : 


/3f\ 

Ujj 


_ I (h -h ) 
h j+1 + h. 2 ch j+l V 



(3.2) 



= 2 


Jbi. 


h. (h.+h 1 ) 

J 1 1+1 


h. h... 
1 1+1 


J±L 


il j+l (h j +h j+l ) 



(3.3) 


Note that the second term of the right-hand side of Eq, (3.2) and 
(3.3) is the "extra error" introduced by the use of a nonuniform grid. 
In general, however, this term is very small if the grid size varies 
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slowly (Blottner, 1974) (this is the case with 3.1). It can be easily 
shown (Blottner 1974) that a variable grid scheme is equivalent to a 
coordinate stretching method if a relation of the form of Eq. (3.1) is 
used to specify both the grid spacing in the variable grid method and 
the relationship between the coordinates for the stretching method. In 
both cases the derivatives are second order accurate in terms of A£ , 
i . e. , 


and 





+ 0(A£ 2 ) 


(3.4)' 



Y 1 
'VYW 


f . 

3 

h. h. 

1 J+l 


+ 


£ 1+i 

h (h +h ) , 
j+1 j j+1 J 


+ 0(AO 

(3.5) 


In the streamwise and spanwise directions the pseudo-spectral method 

9 3 9 ^ 

is used for the calculation of partial derivatives — , , etc. 

dx d z 9 z 

For a given number of grid points, the maximum accuracy is achieved by 
using this method (see Moin et al.„ 1978, for a discussion of the 
accuracy of numerical differentiation operators in terms of modified wave 
number concept). For periodic boundary conditions, which are of interest 
in x and z directions, we can represent a flow variable such as u 
by a discrete Fourier expansion 


where 


u (Xi,x 2 ,x 3 ) 


EE u(k 15 x 2 ,k 3 ) 

n l n 3 


i(k l x l +k 3 X 3 ) 

e 


(3.6) 


k. 

l 


2tt 

(Nh) . n i 
3 


wave number in the x. direction 

J 

number of mesh points in the j direction 
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N. N. 

X x 

— - 2 2" — 1 

h. = mesh size in the x. _ direction . 

^ ' ’ J 

The sum extends over all and • Suppose we wish to compute 

3u/3x-^ ; we may regard (3.6) as an interpolation formula, treating x^ 
as a continuous variable, and differentiate to obtain 


3u XT' V — i(k-x.+k_x,) 

^ “«VV k >ik e 113 3 (3 .7: 

1 n x n 3 

Multiplying both sides of (3.7) by exp (— ikjx.^ — ik^x^) , summing over 
all x^ and x^ , and using orthogonality, we get: 


- ik x u(k 1 ,x 2 ,k 3 ) 


(3.8) 


Thus, in order to compute 3u/3x^ , we simply have to Fourier 
transform u in the x^-direction, multiplying it by ik^ , and take 
the inverse transform of the result; this is called the "pseudo- 
spectral" approach (Orszag (1972), Fox and Orszag (1973)). The use of 
pseudo-spectral method in x and z directions, partially addresses 
the grid resolution problem in these two directions. 

For a limited number of problems with nonperiodic boundary condi- 
tions we can use some other set of orthogonal functions rather than 
{e^^x} (see Orszag, 1971). For completeness and for later use in this 
report, we conclude this section by describing the numerical differ- 
entiation using Chebyshev polynomials. 

' We can express a variable such as f(y) by a discrete Chebyshev 
expansion 

N •• 

= X) T (y) (3.9). 

n=0 


22 



where T (y) is the nth order Chebyshev polynomial of the first kind, 
and double prime denotes that the first and last terms are taken with 
factor \ . Similarly, we can express the derivative of f , which is a 
polynomial of degree N-l , in terms of T (y) . 'We then write 


9f_ 

3y 


N-l, 

12 b T (y) 
„ n n 


(3.10) 


and seek to compute -the coefficients b in terms of a . It can be 

n n 

easily shown (see Fox and Parker, 1968) that the coefficients b are 

n 

given by the following recurrence relations : 


b n-l b n+l 2n a n nl,2,...,N2 

\-2 = 

b H— 1 ' N Vl (3 ' 11) 


Finally, we note that 


T (cos 0) = cos n6 (3.12) 

n 

Thus, the transformation (y = cos 0) which is roughly adequate for 
boundary layer coordinate stretching, renders the evaluation of the 
Chebyshev expansion coefficients, , particularly simple with the 
use of FFT routines. 


3.3 Fundamental Numerical Problem 

In this section we describe an inherent numerical problem associ- 
ated with the fully explicit solution of the dynamical equations in 
primitive form in a bounded domain. Consider the momentum equations 



3t 



+ H. 


(3.13) 
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where IL contains the viscous and convective terms. In the fully 
explicit (time advancing) numerical solution of (3.13) one normally 
specifies an arbitrary initial solenoidal velocity field satisfying the 
no-slip- condition. * Then, one proceeds to solve the appropriate 
Poisson equation for pressure obtained from the application of the 
divergence operator to the momentum equations to ensure that 
V*u = 0 . The resulting pressure is then used together with the com- 
puted H. in (3.13) to advance u^ in time. The Neumann boundary 
condition. 


= vn * V 2 u (3.14) 

dn 

*v 

is normally used in conjunction with the Poisson equation for pressure. 
Here n is a unit vector normal to the solid boundary. This condition 
is obtained from the normal momentum equation evaluated at the solid 
boundary. 

With regard to the boundary treatment, one has two choices: 

a) Enforce the no-slip condition, and time advance the velocity 
field via (3.13) only in the interior domain (not at the boundaries); 

b) Time advance the velocity field throughout (interior domain as 
well as boundaries) . 

If one chooses (a) ; for the tangential momentum equations to be 
satisfied at the boundaries, the initial field would have to be such 
that the p it generates satisfies the Dirichelet condition 

|£ = vx • V 2 u (3.15) 

9t 

(t is a unit vector tangent to the solid boundary) . The momentum equa- 

*v 

tions in the directions tangential to the solid boundary will not 
necessarily be satisfied if the only constraints on the initial field 
are that it be solenoidal and satisfy the no-slip condition. Since the 
tangential momentum equations are not in general satisfied at the solid 
boundary, the Poisson equation will not be satisfied there either, and 
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hence we conclude that in case (a) the continuity equation will not be 

g 

satisfied at the boundary, ( — (n*u) ^ 0) . This can cause serious 

dn ~ ~ 

numerical instability. 

On the other hand, if one chooses case (b) , continuity will be 
satisfied everywhere, but the no-slip condition may not be satisfied, 
and this is unacceptable. 

It should be noted that, if one uses the Dirichlet condition (3.15) 
as the pressure boundary condition then the Neumann condition (3.14) 
will not necessarily be satisfied and hence similar problems will arise 
in either approach (a) or (b) . 

In Appendix B we formally demonstrate the numerical problems 
addressed in this section. In addition, in Section 3.6 it will be 
shown that the numerical problems discussed here can be avoided if one 
uses three-point finite differences to approximate partial derivatives 
in the direction normal to the boundaries. 

3.4 Consistency Conditions for the Initial Velocity Field 

In this section, we present a set of consistency conditions for 
the initial velocity field of the channel flow such that the Neumann and 
Dirichlet problems for the pressure have the same solution, i.e., we 
solve the problem addressed in Section 3.3. 

Fourier transforming the Poisson equation in the streamwise and 
spanwise directions, we get: 


2 ^ 

d P 2~ 

- k z P = 

dy 


(3.16) 


The consistency condition requirements conflict with the proven 
existence and uniqueness theorems for the Navier-Stokes equations. 
Therefore, we emphasize that the problems addressed in the previous sec- 
tion are purely numerical and mathematically there is no difficulty. 
Saffman (P. G. Saffman, 1978, private communication) points out that 
the fact that the Neumann problem does not satisfy the Dirichlet condi- 
tion appears in the nonanalyticity of V^u on the boundary at t = 0 , 
which can be interpreted physically as an~ initial vortex sheet diffusing 
from the boundary. 
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2 2 2 

where k = k^ + , and and k^ are the wave numbers in 

streamwise and spanwise directions respectively. Here, 




8x. 3x. U i U j 
i j 


For k ^ 0 , the general solution of (3.16) is: 


P = cj>(y) + c sinh ky + cosh ky 


(3.17) 


where: 


•Ky) 


r j * i 

r j ~ i 

/ q cosh kn dn 

sinh ky - / Sjishkn ^ 

A 

L4 -1 


cosh ky 


and, c^ and c^ are constants. Thus, for the Dirichlet and Neumann 
problems, we can determine c^ and separately to get P^ and P^ 

which are the solutions of Dirichlet and Neumann problems respectively. 
Note that for the Dirichlet problem to have a solution, we must have 


2 

3 P 
8x3 z 


y=±l 


3 2 P 

3z8x 


y-±l 


The above condition is equivalent to n*V u = 0 on the boundaries 
(y = ±1) , or 


or 


3 

3 2 u 


_3_ 

C5 2 
a W 


3z 

3y 2 

y=±l_ 

9x 

2 

3y 

y=±l_ 


- (3.18) 


ik 3 ^(±1) = i^ h 3 (±d 
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where 


3 2 u 3 2 w 

H 1 = V 72 and H 3 = v Ti- 
dy 3y 

* 

and a) is the vorticity vector. 

Equating and P^ (after some algebra) we arrive at the follow- 

ing constraints for the initial velocity field: 


H 3 (l) - h 3 (-i> 
ikl 


h 3 (i) + h 3 (-i> 


ik„ 


- <KD 

- <KD 


tanhk 

k 

cothk 

k 


H 2 (l) + H 2 (-l) - <fr*(l)J 
H 2 (1) - H 2 (-1) - 4>* CD^] 


(3.19) 

(3.20) 


Therefore, for a successful, fully explicit numerical simulation, the 
initial velocity field must satisfy the following conditions: 

• it must be solenoidal, 

• it must satisfy the no-slip condition, and 

• it must satisfy (3.18), (3.19), and (3.20). 

Note that for k 3 = 0 and k^ ^ 0, one can use (3.19) and (3.20) with 
the subscript 3 replaced by 1. 


3.5 Conservation Properties 

As was pointed out by Phillips (1959) , numerical integration of the 
finite-difference analog of the Navier-Stokes equations may introduce 
nonlinear instabilities if proper care is not taken. Differencing the 
transport terms in the form of (2.5) will automatically conserve momentum 
in an inviscid flow. However, in general, the computation becomes un- 
stable and the kinetic energy increases. This can happen in spite of the 
dissipative nature of t. . and the viscous terms. The nonlinear insta- 
bility arises because the momentum conservative form does not necessarily 
guarantee energy conservation (in the absence of dissipation), and the 
effect of truncation errors on the energy is not negligible. 

Moin et al. (1978) have shown that writing the dynamical equations 
in the form of (2.6) results in vorticity, momentum, and energy conser- 
vation for a large class of differencing schemes. Therefore, in all the 
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calculations reported here, we use the dynamical equations in the form 
shown by Eqn. (2.6). 


' 3;6 ' Explicit Time Advancing' 

By introducing one plane of grid points just outside of each bound- 
ary, one is able to obtain some degree of freedom. With proper use of 
this freedom, one can avoid the problem discussed in Section 3.3 (case a) 
The reader should be cautioned that here we are strictly referring to the 
explicit numerical solutions in which three point finite differences are 
used for the numerical differentiation. (However, the latter statement 
does not apply, for example, to the cases in which Chebyshev polynomials 
are used in a finite series expansion to represent a flow variable and 
its derivatives in the normal direction (see Sec. 3.2).) In practice, 
one can determine the normal velocity at the exterior point such that 
the continuity equation evaluated at the wall. 


ov 

3y 


= 0 

y=±l 


(3.21) 


is identically enforced. This velocity, in turn, is used in obtaining 

the Neumann boundary condition for pressure. For the proper choice of 
2 

the numerical V operator for the Poisson equation, the reader is re- 
ferred to Moin et al. (1978). 

For explicit time advancement, a second-order Adams-Bashforth method 
was used. It has been shown by Lilly (1965) that this method is weakly 
unstable, but the total spurious computational production of kinetic 
energy is small. The Adams-Bashforth formula for u^ at time step 
n 4- 1 is 

^J +1 = uj + At (f^ - j-Cj"* 1 ) + 0(At 3 ) (3.22) 


where 




/3u. 3u.\ 3 t.. . 3 2 u. 

— I a _ j \ _ 3P _ ii <. _1 _x_ 

U j \3x. 3x. / 3x. 3x. il Re 3x.3x. 

v 1 1/ i 1 T J 3 
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Using the above method, we have successfully integrated the governing 
equations for the numerical simulation of turbulent channel flow (not 
reported here). However,^ due to the presence of a very fine mesh near 
the boundaries, one is forced to use extremely small time steps. This 
stringent requirement is caused by the well-known numerical stability 
criterion of the diffusion equation. 


3.7 A Semi-Implicit Numerical Scheme 


As was mentioned in the previous section, due to the presence of 
diffusion terms in the governing equations, the time-step requirement of 
a fully explicit method becomes severe. To circumvent this difficulty, 
we have devised a semi-implicit algorithm. All the results reported here 
were obtained using this method. Thus, in what follows, we outline a 
method which treats part of the diffusion terms and pressure in the dy- 
namical equations implxcitly, and the remaining terms explicitly. The 
equation of continuity is solved directly. 


Let us start with Eqn. (2.19), written in the following form: 


where 


H. 

l 


9u . 

i 

9t 


H. 

l 


9P 

9x. 

i 


(Vt + rH 


Q 2~ 
9 u. 


9 x. 


(no summation) 

(3.23) 



<k 



(no summation) 


C. 

i 


1 + 6 


i2 
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d. = 1 - S.„ 

i x2 


V T ■. < y*r IE 2 , ?3. > - > x 1 ,*3 


< > indicates the average of bracketed quantity in x 1 -x^ plane, 

xi, x 3 1 J 

v T = v m - v m 

T T -T 

The rationale for this decomposition of will be explained later in 

this section. For time advancing, we are going to use the Adams- 

Bashforth method (see Sect. 3.6) on H. , and the Crank-Nicolson method 

1 ' 2 — 2 

(Richtmyer and Morton, 1967) on 9P/9x.. and 9 > in the right- 

hand side of Eqn. (3.23). For convenience, we evaluated V T at time 
step n. Thus, we have: 


—n+1 — n . / 3 n 1 TT n-l\ At /3P n+1 3P n \ 

u i = u i + At [j H i " 2 H i ) ~ T [~^r + a* l } 


+ c.v n ) 

\Re x i T/ 


( r. 2— n+1 

3 u. 


8x 2 


B i (x 2 ) = - 


2/At . 

(sr* 0 ! 3 ? 6 ^) ’ 


rearrangement of Eqn. (3.24) yields: 

cs 2— n+1 , - 

i „ -n+1 „ At 3F n+ 

YY + S i u i * S iY — 

dX 0 dX. 

I 1 


■ e i“i + B i At (f H 1 " I H i _1 ) 


~2— n 

, n d u . 

a At 9P 1 

" P i 2 9x. “ Q 2 
1 dx^ 


(no summation) 


Finally, we write the continuity equation at time step n+1: 


1 


= 0 


(3.26) 
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Now let us Fourier transform Eqns. (3.25) and (3.26) in x-^ and 
Xg directions. This transformation converts the set of partial differ- 
ential equations (3.25) and (3.26) to a set of ordinary equations for 
every pair of Fourier modes k^, with x 2 as the independent vari- 
able. Note that the dependent variables have superscript n + l. In the 
remainder of this section all the dependent variables are to be inter- 
preted as two-dimensional Fourier transformed quantities. Fourier trans- 
forming equations (3.25) and (3.26) results in the following set of 
ordinary differential equations for the dependent variables: 


3 2 u n+1 

1 , n n+l . ^ At „n+l 

~rr- + h u i + lk A x p 

3x_ 


<>2 n+l 
3 u^ 

3x? 


4. a ,, n+1 + R ^ 3P 

+ B 2 u 2 + 3 2 2 3x. 


n+l 


~2 n+l 
3 u 3 

3*5 


, c n+l . „ At n+l 

+ e 3 u 3 + ik 3 B 3 — p 


r. n+l 

•V «+l 3U 2 

lk l"l + 

and k„ 


h u ? + 6 


- ik 


- 6 2 X 


- ik 3 B 3 


+ ik 3 u 3 


n+l 


, At J 
*1 2 1 

3H? - H?' 
^ 1 1 

At 

„2 n 
3 u 3 

2 P 

3x 2 

CM 

( 3H 2 - H 2 

3p n 

3 2 ," 

3x 2 

3x 2 

s At , 
3 2 1 

( 3H 3 - H"' 

At 

r.2 n 
3 u 3 

2 P 

9x 2 

■ 0 



(3.27a) 


(3.27b) 


(3.27c) 


(3. 27d) 

we have four coupled linear ordi- 


Thus, for every pair of k^ ...... 

nary differential equations with u^ + 1 (k 1 ,x 2 ,k^) , u 2 +1 (k^ ^ 2 ,^) , 

U 3 + 1 (ki’X 2 >k 3 ), and P n+ ^ (k^ ,x 2 ,k 2 ) as unknowns. Note that, with no 
further^ complications, one can treat more terms in Eqn. (2.19) (e.g., 

1 3 1 3 2 ...... 

—7 , etc.) implicitly. 


Re 


3x 


Re 
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Finally, it should be mentioned that, in order to avoid evaluating 
complicated convolution sums, we have decomposed to its planar 

average, V (y) and "fluctuating" component V* (x.. ,x 9 ,x~) . We have used 

i 2 2 x ~ J 2 , 2 . 

expii'ci't time advancing for V.J, (3 u^/3x 2 ), w ^ ereas V T (3 u -j/^ x 2 ) 
advanced by a partial implicit scheme. This decomposition of V T may 
not be an optimum one from the standpoint of numerical stability and 
accuracy. Other -choices are possible. For example, one can decompose 
V as follows : 

V t (x 1 ,x 2 ,x 3 ) = max (v T ) + V^x-^x^x^ 

X 1’ X 2 >X 3 

Although we did not incorporate any other decomposition than the one used 
here, relatively simple numerical experiments with the diffusion equation 
may result in a better decomposition for V^. 


3.8 Finite-Difference Formulation and Boundary Conditions 


In order to solve Eqns. (3.27) numerically, we use the finite dif- 

2 2 

ference operators (3.2) and (3.3) to approximate 3/3x 2 and 3 /3x 2 . 
Having done this, we shall have a set of linear algebraic equations for 
the Fourier transform of the dependent variables. This system of alge- 
braic equations is of block tri-diagonal form and can be solved very 
efficiently. However, in order to close the system we must provide a 
set of boundary conditions, i.e., we have to specify the values of u^, 
u 2’ U 3 J an ^ ^ at t * ie so -^-^ boundaries. 

Implementation of velocity boundary conditions poses no problem; we 
simply set the value of the velocity vector at zero on the walls. In 
order to obtain the pressure boundary conditions, we note that evaluation 
of Eqn. (3.27b) at the solid boundaries yields: 


*2 n+1 
3 u. 


3x 2 


+ 3 


At 3P n+1 
2 2 3x 2 


x 2 =±l 



At 3P^ 
2 2 3x 2 
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.2 n 
3 u 2 

a 2 

3x„ 


x 2 ±;L 


Consider the following Neumann boundary condition for pressure: 



x 2 =±l 


-2 

1 3 U 2 


Re 


T' 8x: 


x 2 =±l 


(3.28) 
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Equation (3.28) was obtained from the Fourier transform of Eqn. (2.19, 
i = 2), and evaluated at the solid boundaries. It is clear that this 
equation is consistent with the numerical analog of that equation (3.27b) 
evaluated at the walls. Note that 


o At 

B 2 T 


x 2 =±1 


- Re 

T 


Thus, we formally use Eqn. (3.28) as the pressure boundary condition. 
However, for closure the finite-difference equations require the value 
of pressure at the boundaries, not its normal derivative. For this we 
use the following difference relation in conjunction with the difference 
analog of Eqn. (3.28): 


where h . 

3 

wall. 


1 / 9P 


2 l9x 0 


+ 


9p 


j-2 3X 2 


P - P 

Ji+k... 1 


j=3j 


h. 

3+1 


3=2 


+ OCh*) 


(3.29) 


x 2 “ x 2 ’ 3 “ 2 indicates the grid point on the lower 

3 3-1 


Substituting the finite-difference analog of Eqn. (3.28) into the 
left-hand side of Eqn. (3.29) and using the finite-difference form of 
the continuity equation at the wall, we obtain: 

30) 


2P„ 


P 2 " 


(h^ + h^) 


2u„ — i 

3 

Re h„h„ 

X 2 3 



+0 


h 4 + “3 


(3. 


An analogous relation is used for the value of the pressure at the 
upper wall (j = N - 1). Note that the pressure is still indeterminate 
by a constant, as it should be due to the use of Neumann boundary condi- 
tions; i.e. , we are not using Dirichlet boundary conditions. 

In the case = 0, a special solution technique must be un- 

dertaken. First observe that in this special case Eqns. (3.27a) and 
(3.27c) are independent of each other and Eqns. (3.27b) and (3.27d). 
Furthermore, the former two equations are of simple tridiagonal form and 

n4-1 

can be solved directly to yield u^ (OjX^^) and u^ (OjX^^). Second, 
the continuity equation together with the boundary conditions for u^ 
yield 
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u 2 (0,x 2 ,0) 


0 


(3.31) 


Since pressure is indeterminate by a constant, let 


P(0,x 2 ,0) 



0 


(3.32) 


Using Eqns. (3.30), (3.31), and (3.32) in conjunction with the finite 

difference analog of Eqn. (3.27b) allows one to solve for P (0,x 9 ,0), 

z j 

j - 3,4,.,., N-KL . J 

Before concluding this section, we emphasize that, in obtaining the 
pressure boundary conditions, we used a momentum equation evaluated at 
the boundary. We were able to do this because the finite difference 
equations are generally enforced inside the spatial domain and not on its 
boundaries. Consequently, we did not use a redundant equation. Consider 
for a moment a hypothetical case in which we have the means to integrate 
the governing equations of motion analytically. In this case, the equa- 
tions of motion are and should be valid at the boundaries as well as in- 
side the domain (we do not have any singularity at the boundaries). So, 
in this case, use of momentum equations for the pressure boundary condi- 
tions will not provide any new information. The roots of this apparent 
dilemma lie in the basic physics of fluid mechanics. The fact is that 
physics does not provide a priori boundary conditions for pressure. 

A manifestation of this dilemma will appear if, for example, Cheby- 
shev polynomials are used in a finite series expansion to represent a 
flow variable in the y direction (see Section 3.2). However, since 
the equation of continuity is solved directly, it appears that the numer- 
ical problems which were addressed in Section 3.2 will not cause any dif- 
ficulty if one uses Chebyshev polynomials in conjunction with the semi- 
implicit scheme developed here. 


3.9 Computational Details 

The numerical solution of the equations described here (see also the 
next chapter) were carried out on the CDC 7600 computer at NASA-Ames Re- 
search Center. The dimensionless time step, during most of the 
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calculations, was set at At = 0.001. Throughout the computations re- 
ported here, the values of the following quantities, 

c^(t) = MaxjAt 

and 

c 2 (t) = Max 

did not exceed 0.3 and 0.08, respectively. In addition, the numerical 
stability was checked by a 200-step numerical integration in which the 
value of At = 0.0005 was used. The computer-generated results of this 
run agreed (within two significant figures) with the corresponding numer- 
ical integration in which the value of At = 0.001 was used. Comparison 
was made at the same total time of integration. 

The computer time per time step was approximately 20 seconds (CPU 
time). However, the present computer program is not an optimum one, and 
we believe that at least a 25% savings in computer time can be achieved 
by some modifications of this program. 

Finally, it should be noted that, in the present computation, approx- 
imately 80% of the small-core memory and only 50% of the available large- 
core memory of the CDC 7600 was used. Therefore, a computation with 
twice as many grid points as the present one is possible using the avail- 
able core memory of the CDC 7600. 
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Chapter IV 


INCOMPRESSIBLE TURBULENT CHANNEL FLOW 


4.1 Physical Parameters 

In order to solve Eqns. (2.19), we need to specify Re^, Reynolds 
number based on channel half-width 5 and shear velocity u^_. In the 
present numerical simulation of turbulent channel flow, Re^ = 640.25 
was used. In their experimental investigation of the mechanics of orga- 
nized waves, Hussain and Reynolds (1975) considered a channel flow with 
the same Reynolds number. The mean flow parameters of their experiment 
are listed below. 

Re = 13800 

u 

= 0.0464 

o 

U 

~ = 0.881 
o 

U = 21.9 ft/sec (6.67 m/s ec) 

o 

where Re is the Reynolds number based on channel half-width, S, and 
the centerline velocity, U q ; U is the mean profile average velocity, 
and u^ is the shear velocity. 

4. 2 Initial Condition 

A number of initial velocity fields were explored. With the simple 
sub-grid scale model used, it is important that the initial turbulence 
field be able to continually extract energy from the mean flow in order 
that a statistically steady solution develop. For this purpose, we em- 
ployed the governing equations of small disturbances used in hydrodynamic 
stability theory (other choices are possible) to obtain a velocity field 
with negative Reynolds stress. 

The equations for a small wave disturbence u^ on a parallel mean 
flow U(y) are (Lin, 1955, Eqn. (1.3.9)): 
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^ . i(ax+Bz-act) , 

u. = u. (y) e + coni 

11 


iau.. + iBu_ + Du 0 = 0 


iwu^ 4- Uiau^ + DU • U2 - - iotP + (D^ - k 2 ) 


icou^ + Uiau^ = - DP + — ^ (D 2 - k 2 ) u 2 


A y\ a ■! 2 2 A 

iuu„ + Uiau„ = - iBP + ^~ (D - k ) u- 

3 3 Re 3 


Here U) = - ac is the (complex) frequency, and D = d/dy. 
The Squire transformation (Lin, Eqn. (3.1)), 


.2 2 q 2 

k = a + 3 


(4.1a) 
(4.1b) 
(4.1c) 
(4. Id) 
(4.1e) 


(4.2a) 


v = a 2 (4.2b) 

au^ + Bu 3 = ku (4.2c) 

A 

permits reduction to a single fourth-order equation for v, the Orr- 
Sommerfield equation (Lin, Eqn. (1.3.15)): 

v (D 2 - k 2 ) 2 v - ictRej(U - c) (D 2 - k 2 ) v - D 2 u • vj (4.3) 

For a given set of a. Re, 3, and U(y), (4.3) is solved numerically 

using the algorithm of Lee and Reynolds (1967). 

After final calculation of v, u is calculated from (4.1b), and 
p is calculated from (4.1c) and (4.1e). The results are then used to 
solve for u^, via (4.1c). Solution of (4.1c) is carried out numeri- 
cally using a second-order algorithm. Starting at the centerline of the 
channel, two solutions, each satisfying the centerline boundary condi- 
tions (here we are primarily concerned with symmetric u 2 and anti- 
symmetric u^ and u^) are constructed using the Kaplan filtering tech- 
nique to maintain linear independence. These two solutions are then 
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combined to satisfy the wall boundary conditions. The eigenvalues are 
automatically adjusted until an eigensolution is obtained. 

For the Reynolds number under consideration (Re^ = 640.25) and 
with proper choice of a, g, and U(y), one can ob.tain a set of u^,- 
u 2 j and u^ such that the corresponding Reynolds stress has the same 
sign as - Du. This corresponds to an unstable disturbance from the view 
of hydrodynamic stability theory. The resulting three-dimensional dis- 
turbance extracts energy from the mean flow in a continuous fashion. In 
the present study we have used a = 1.0, 3 = 1.5, and the mean velocity 
profile: 


U(y) = 10(1 + cos Try) 

for the generation of initial disturbances. 

This profile was chosen with due consideration to the proper repre- 
sentation of the resulting disturbances on the grid system in the normal 
direction. In addition, note that the above mean velocity profile has 
inflection points (at y = ± y) which produces Kelvin-Eelmholtz type 
instability. 

In order to avoid a net momentum in the spanwise direction, one can 
add two oblique waves with the same amplitude that are traveling in the 
directions which are at angles of <f> and - <p with the streamwise, x, 
direction. Combining two oblique waves in this fashion yields a set of 
streamwise vortices (roll cells). Thus, the following velocity field was 
used as the major part of the initial disturbance (initial large eddies) : 

u^tx^jz) = Afu^Cy) cos gz e 1CCx + conj ] 

u 2 (x,y,z) = A[u 2 (y) cos gz e 10tx + conj] 

u^x^z) = Aju^y) sin gz e iax + conj] 

A. 

Here, A is a constant, a = 1.0, g = 1.5, and u^(y) (i = 1,2,3) are 
the eigensolutions of the linearized equations. In order to allow the 
development of all the waves that can be resolved on the grid system, a 
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solenoidal velocity field with random phase was added to the above veloc- 
ity field. Furthermore, to ensure the initial dominance of the u. 

i 

field, the amplitude of random field was about 10% of the maximum apli- 
tude of u^. Finally, in order to avoid a very long time numerical in- 
tegration, the measured mean velocity profile of Hussain and Reynolds 
(1975) was used as the initial mean velocity. 

4 . 3 Preliminary Numerical Experiments 

In the following three sections we shall present and discuss various 

calculated quantities pertinent to turbulent channel flow. The results 

will consist of running time averaged mean velocity profile and turbulence 

statistics, horizontally (xz plane) averaged turbulent quantities, and 

some instantaneous velocity profiles. However, first, it is instructive 

to discuss some of our initial numerical experiments (failures). 

In our first integration attempt, we observed that the absolute 

value of the horizontally averaged Reynolds stress, < uv > , decreased 

continuously in time. This vanishing trend occurred in spite of the fact 

that the Reynolds stress profile was below the expected value. The total 

time of integration was approximately 1 nondimensional unit, and the 

value of eddy viscosity constant, C g , was specified to be 0.2 (see Moin 

et al., 1978). It is interesting to note that the profiles of 

<(u~<u>) > were generally increasing, and the corresponding 

—2 1/2 

profiles of < v > were decreasing slightly . In other words, the 

correlation between (u - < u > ) and v, and not the respec- 
tive intensities, had a rapid vanishing trend. At this point it was de- 
termined that the effective Reynolds number (taking the eddy viscosity 
into account) was probably too small for a small amplitude disturbance to 

grow. With this in mind, and noting that the production of Reynolds 

—2 

stress is directly proportional to < v > , the existing turbulent vel- 
ocities were multiplied by a factor of two (and the Reynolds stress was 
amplified by a factor of four). Note that no changes were made to the 
final mean velocity profile, < u > . In fact, at this time < u > was 
deviated considerably from its original profile. 


39 



Using the resulting velocity field as a new initial condition (in. 
what follows, we shall call this velocity field "field A"), we carried 
out two parallel computations, one with C s = 0.44 and the other with 
C s - 0.2. In the former case, the Reynolds stress profile grew contin- 
uously for a nondimensional time, t, of 0.3. However, during a further 
integration period (t = 0.7), it decayed drastically to a vanishing 
level. Thus, it was concluded that the value of 0.44 for the sub grid 
scale model constant is too large, causing turbulent motions to damp out. 

The results to be presented in the following sections were obtained 

using the value of 0.2 for C . This value is probably not the omptimum 

s 

one (more likely the optimum value is between 0.2 and 0.3); however, in 
the absence of a more rigorous sub grid scale model formulation, further 
adjustments of C g seem to be unjustified. 

4.4 A Time History of the Horizontally Averaged Turbulent Quantities 

As was pointed out in the previous section, we use the velocity field 
A as the new initial condition. Fig. 4.1 shows the horizontally averaged 
resolvable shear stress < uv > of this field. For purposes of discus- 
sion, we concentrate on the lower half of the channel in this section. 
Furthermore, due to the relationship between the materials to be discussed 
herein and the bursting process in a turbulent boundary layer, virtually all 
of our discussion will be concerned with the region near the (lower) wall. 

Figure 4.2 shows the < uv > profile at the non-dimensional time , 
b = 0.45. It can be seen that the resolvable shear stress profile has 
increased considerably. In particular, near the wall it has increased 
significantly beyond the expected equilibrium (time- averaged) value. Figs. 
4.3, 4.4, and 4.5 show the profiles of the same quantity ( < uv > ) at 
three later times (t = .65, .85, 1.05, respectively). It is clear that, 
especially in the region -.95 < y < -.7, a dynamic process exists which 
nearly repeats itself in time. If we carry out the integration still 
further, we see the same behavior (almost cyclic) in the < uv > profile. 


* 

One nondimensional time unit corresponds approximately to the time 
in which a particle moving with centerline velocity travels 22 channel 
half -widths. 
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Figs. 4.6 and 4.7 show the vertical distribution of < uv > obtained 
at two later times corresponding to t = 1.425 and t = 2.025, respec- 
tively . 

Since the production of the resolvable turbulent kinetic energy is 

directly proportional to < uv > , it should be interesting to study the 

— — 2 1/2 

effect of the cyclic behavior of < uv > on < (u - < u > ) > 

— — 2 1/2 

Figs. 4.8, 4.9, and 4.10 show the profiles of < (u - < u > ) > in 

+ — 
the vicinity of the wall (y < 128) . They correspond to the < uv > 

profiles presented in Figs. 4.5, 4.6, and 4.7, respectively. Examination 

of these figures shows clearly the effect of production on the 

— — 2 1/2 

<(u-<u>) > profile. It can be seen that, during the times at 

which < uv > has a relatively high value, the corresponding 

— — 2 1/2 

< (u - < u > ) > profile possesses a pronounced local maximum. It 

is interesting to note that, during the quiescent (low < uv > ) periods, 
the turbulence energy level is still quite large. In fact, a close exami- 
nation of Figs. 4.9 and 4.10 reveals that, during these times, the energy 

that gave rise to the local maxima is distributed throughout the 
2 1/2 

< (u — < u > ) > profile. This results in a wide maximum (in con- 

— — 2 1/2 

trast to a sharp local one) in<(u-<u>) > 

During their investigation of the ''bursting" process in a turbulent 
boundary layer, Kim et al. (1968) showed that, while the bursting process 
indeed contributes to' the turbulent energy, its main effect is to provide 
turbulence with u' and v' in proper phase to give the large turbulence 
stress required for an increase in production. This is precisely what is 

observed here. To clarify this point, consider, for example, Figs. 4.6 

q. 

and 4.7. If we focus out attention on the vicinity of y =64 (y = -.90), 
we see that the value of < uv > in Fig. 4.6 is about twice the corres- 
ponding value in Fig. 4.7. On the other hand, the corresponding value 
— — 2 1/2 

of < (u - < u > ) > in Fig. 4.9 is only 6% higher than the one in 

—2 1/2 

Fig. 4.10. And the corresponding values of < v > (Fig. 4.11) and 
2 1/2 

< w > (Fig. 4.12) show no significant change during this period. 

—2 1/2 

This is expected, since the governing equations of < v > and 
2 1/2 

< w > do not contain direct production terms. These quantities can 
only be fed by the inter-component transfer mechanism, which is generally 
a slow process. 
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We conclude this section by considering, once again, our initial 
numerical experiment (see Section 4.3). Recall that, during the first 
integration attempt, < uv > had a rapid vanishing trend while the in- 
dividual components < (u - < u. > ) >- -and- -< v > -did" not" (the " ' 

latter had a slight decreasing trend) . With this and the discussion of 
the present section in mind, one can see the importance of the phase re- 
lationship between (u - < u > ) and v. Indeed, the correlation between 
(u - < u > ) and v is the essential factor for the maintenance of tur- 
bulence. We believe (on the basis of a cursory scan) that the increase in 
< uv > is also highly localized in space. 

It should be noted that, in a computation with a large number of 
mesh points in the horizontal planes, the transitory behavior of < uv > 
described in this section, will not occur. In this case, the horizontal 
averaging is approximately equivalent to long-time averaging; and in order 
to study the relationship of the bursting process to the turbulence stress, 
one should study the time history of the (u - < u > ) v profile at one 
(x,z) location. Such a study, in turn, would yield the-mean bursting 
frequency. 

4.5 Detailed Flow Structures 

In this section we examine some of the detailed flow patterns. Par- 
ticular attention will be given to instantaneous velocity profiles. Fig. 
4.13 shows typical instantaneous streamwise velocity profiles, u. These 
profiels are obtained at the same location (x = 0, z = 13 h^), but at 
two different times (t = 1.625, t = 1.825). For comparison, the mean 
velocity profile is also included. Fig. 4.14 shows the corresponding 
normal velocity profiles, obtained at the same location and times. Exam- 
ination of these figures reveals that the profile with a momentum defect 
(with respect to the mean) corresponds to a case in which fluid is being 
ejected from the wall (v > 0), while the profile with excess momentum 
corresponds to a case where the flow is toward the wall (v < 0). In 
addition, both pairs (((u - < u >) > 0, v < 0) and ((u - < u >) < 0, 
v > 0)) have positive contributions to the resolvable Reynolds stress 
and, hence, they contribute to the production of turbulence. 
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The velocity profiles presented here are in good qualitative agree- 
ment with the flow visualization data of Kim et al. (1968) and Grass (1971). 
In their study of the bursting process in a turbulent boundary layer, Kim 
et al. observed that during the gradual lift-up of low speed streaks from 
the sublayer, inflectional instantaneous velocity profiles were formed. 

In fact, the appearance of the inflectional profile was used as one of 
their criteria for the detection of the bursts. 

Using the terminology of Grass, the u profile with momentum defect 
corresponds to the ejection phase of the bursting process while the profile 
with excess momentum corresponds to the inrush phase (sweep). In the 
lower left-hand corner of Fig. 4.13, we have included the ihstantaneous 
velocity profiles from the measurements of Grass (1971) in a flow over a 
smooth flat plate. In Figs. 4.15 and 4.16, the same quantities as in Figs. 
4.13 and 4.14 are plotted, but they are obtained at a different location 
and at different times (x = 10 h^, z = 10 h^> t = 1.05, 1.275). The 
same behavior (qualitatively) as in Figs. 4.13 and 4.14 are displayed by 
Figs. 4.15 and 4.16. Fig. 4.17 shows the instantaneous streamwise veloc- 
ity profiles obtained at time t = 2.025, but at two different (x,z) 
locations. This figure, together with Figs. 4.13 and 4.15, clearly demon- 
strate the highly three-dimensional and unsteady nature of this flow. 

The reader is cautioned against establishing a direct relationship 
between the times, t, at which the instantaneous profiles are presented 
here, and the corresponding times at which < uv > assumes a relatively 
high or low value (see the previous section) . Recall that in this section 
instantaneous velocity profiles were presented at one (x,z) location, 
while in the previous section we were concerned with the planar averages 
of < uv > . At most we can say that, during the times at which < uv > 
has a relatively high value near the wall, there are more locations where 
the relationship between the u and v profiles are the same as those 
shown in Figs. 4.13 and 4.14 (((u - < u >) > 0, v < 0) or ( (u - < u >) 

< 0, v > 0)). This is in contrast to the times at which < uv > has a 
relatively low value. 

At this point, let us consider the spanwise instantaneous velocity 
profiles. Figs. 4.18 and 4.19 show a typical spanwise variation of the 
streamwise velocity u in the vicinity of the lower wall (second grid 
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point away from the wall, y 53 3.85) at eight consecutive strearawise 
locations. The profiles presented here are obtained at time t = 1.05. 
These figures demonstrate distinct regions of high-speed fluid located 
adjacent to the low-speed ones. In addition, these profiles clearly show 
the long strearawise extent of the high- and low-speed streaks. In. their 
visual studies, Runstadler et al. (1959, 1963) (see Section 1.2) demon- 
strated that the viscous sublayer consists of relatively coherent struc- 
tures of low- and high-speed streaks alternating in the spanwise direction 
over the entire wall. '‘It appears, therefore, that at least there is a 
qualitative agreement between the calculated results and the laboratory 
observations. Figs. 4.20 and 4.21 show the spanwise profiles of u at 
the same locations as in Figs. 4.18 and 4.19, but at time t = 1.425. 

Once again, these profiles show the coherent structures of alternating 
low- and high-speed streaks. Note that the profiles shown in Figs. 4.20 
and 4.21 are generally different in magnitude and details of structures 
from those presented in Figs. 4.18 and 4.19 (see, for example, the pro- 
files at (x = 0, and x - 4 h^). Fig. 4.22 shows typical spanwise vari- 
ation of v and w, obtained at y + = 3.85, t = 1.05, and x = 4 h^. 
The rapid spanwise variations of v and w clearly show the lack of grid 
resolution in the z direction (see the following discussion) . Neverthe- 
less, these profiles demonstrate, once again, that the viscous sublayer is 
the region of high flow activity, and it is three-dimensional. In addi- 
tion, the spanwise variations of v indicate the distinct presence of 
secondary longitudinal vortices in the wall region. 

Before concluding our present discussion of the spanwise velocity 
profiles , it is appropriate to make a comment about the grid resolution. 
Examination of the spanwise velocity profiles, in particular v and w, 
seems to show that a better resolution in the z direction is required 
(see Section 3.1 and also note that our streak spacings are far larger 
than experimental observations). In other words, more grid points in the 
spanwise direction are necessary to represent the relatively rapid varia- 
tions of the velocities (streaks) properly. This is necessary in spite 
of the fact that the pseudo -spectral method is used for numerical differ- 
entiation in the z direction. However, since the eddies away from the 
boundaries are larger than the ones near the walls (see Fig. 4.23), 
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it is probably sufficient to have more grid points just in the vicinity 
of the walls. This requires a non-re ctangular grid system (conical), 
which is generally accompanied by computational difficulties. Finally, 
Fig. 4.24 shows typical streamwise variations of v and w which are 
obtained at t = 1.05, y + = 3.85, and z = 8h^. Note that, in spite 
of the fact that these profiles are obtained at the same plane as those 
in Fig. 4.22, the streamwise grid resolution seems to be adequate. How- 
ever, it appears that the streamwise extent of the computational box, 

L , is too small, 
x 


Running Time Average of Mean Velocity Profile and Turbulent Statistics 


In this section, we shall present the calculated mean velocity pro- 
file and turbulence quantities, averaged over horizontal planes and in 
time. The total averaging time is about one dimensionless time unit, which 
is much smaller than corresponding time intervals commonly used in labora- 
tory measurements. However, the horizontal averaging should somewhat im- 
prove the overall statistical sample. In addition, note that, during the 
time interval used for the averaging (1.05 < t < 2.025), the resolvable 
shear stress profile < uv > traversed (roughly) one cycle (see Sect. 4.4). 

Vertical profiles of the resolvable mean Reynolds stress, < uv > 
(unless otherwise stated in this section, < > indicates horizontal as 
well as time averaging) , and the total Reynolds stress 


< uv > + < - V T 



> 


are shown in Figs. 4.25 and 4.26a. These profiles indicate that an approx- 
imately steady mean velocity is obtained. In other words, the average 
Reynolds stress profile has nearly attained the equilibrium shape which 
balances the downstream mean pressure gradient in the regions away from 
the walls. In the vicinity of the walls, the viscous stresses are signi- 
ficant, and they, together with the total Reynolds stress, balance the 
mean pressure gradient. Moreover, it should be noted that the subgrid 
scale cobtribution to the total Reynolds stress is significant only in 
the vicinity of the walls (see Figs. 4.25, 4.26a, and 4.26b). 
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Figure 4.27 shows the profile of < u >, the mean velocity, aver- 
aged over both halves of the channel. The latter averaging was performed 

* 

in order to improve the overall statistical sample . The calculated mean 
velocity profile shows a distinct logar-ithmic region. In addition, the 
agreement with experimental data is satisfactory. 

Figures 4.28, 4.29, and 4.30 show the profiles of the resolvable and 

kk 

total turbulent intensities averaged over both halves of the channel 
The contribution of the subgrid scale motions to the turbulent intensi- 
ties is obtained from Eqn. (2.8) and from 


~ < u'.u ! > = 

3 ix 


< f v 2 /(C £) 2 > 


(4.4) 


C = 


.094 


(see Moin et al., 1978, or Lilly, 1967). 

It should be noted, however, that due to the presence of a rela- 
tively coarse grid and the high degree of anisotropy in the channel flow, 
the validity of Eqn. (4.4) is questionable, especially in the vicinity of 
the walls. For comparison, we have also included some of the available 

experimental data in Figs. 4.28, 4.29, and 4.30. Examination of these 

—2 1/2 

figures reveals that, asidejfrom a relatively high value of < v > 
in the vicinity of the channel centerline, the qualitative behavior and 
the relative magnitudes of the turbulent intensity profiles are in accord 
with the experimental measurements. The quantitative agreement of calcu- 
lated turbulent intensities with experimental measurements is good for 

2 1/2 ' 2 1/2 
and < w > and fair for < v > • 


. , . s> 2 1/2 
< (u - < u >) > 


•One may note that the subgrid scale contribution to the total stream- 

wise and spanwise turbulent intensities is relatively small. However, 

Fig. 4.30 shows that, especially in the vicinity of the walls, a large 

2 1/2 

fraction of the vertical turbulent intensity component < v > lies 


The maximum deviation of the calculated mean velocity profile in each 
half of the channel from the one presented in Fig. 4.27 is less than 5%. 

'k'k 

The maximum deviation of the calculated turbulent intensities in each 
half of the channel from the ones shown in Figs. 4.28, 4.29, and 4.30 is 
less than 12%. 
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in the subgrid scale motions. The deficiency in the contribution of the 

—2 1/2 

resolvable motions to < v > suggests that a sub grid scale model 

—2 1/2 

which extracts less energy from < v > might be required. This, in 
turn, may necessitate the use of transport equations for the subgrid 
scale Reynolds stresses (Deardorff, 1973). 

For many problems in fluid mechanics, a knowledge of pressure fluc- 
tuations is desired. For instance, the generation of noise by turbulence 
is related to the distribution of pressure fluctuations. In addition, 
information about the structure of turbulence in the vicinity of the wall 
may be gained from the knowledge of pressure fluctuations at the wall. 
Unfortunately, due to experimental difficulties, direct measurements of 
pressure fluctuations within a turbulent flow are not possible. However, 
from experimental measurements and theoretical considerations, a number 
of investigators have obtained values for the root-mean-square wall pres- 
sure fluctuations in a turbulent boundary layer (see Willmarth and Wool- 
ridge, 1962, and Lilley, 1960). 

In our computer runs, we neglected to calculate the running time 

average of the RMS wall pressure fluctuations. However, we had stored 

the pressure and velocity fields at several dimensionless times. Table 

4.1 shows a time history of root-mean square value of the resolvable wall 

—2 1/2 

pressure fluctuations, < p > /t . Here, < > indicates the aver- 

age of the bracketed quantity over all the grid points on a wall. 


Table 4.1 

RMS Value of Wall Pressure Fluctuations 


Dimens ionless 
Time , t 

< P 2 >in i*v 

Lower Wall 
(y = -1) 

< 7 2 > 1/2 /^ 

Upper Wall 
(y = +1) 

1.05 

2.04 

2.01 

1.275 

1.78 

2.81 

1.425 

1.87 

2.50 

1.625 

2.13 

2.00 

1.825 

2.00 

1.95 

2.025 

1.72 

2.01 
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The average value of the entries in this table (an approximation 

—2 1/2 

for the running time average), < p > /t =* 2.07, is in accordance 
with experimental measurements (see Wilimarth and Wooldridge, 1962, for 
the data from several measurements-) and theoretical estimates (Lilley, 


1960). 

A quantity of particular interest to turbulence modelers is the 

3 

pressure work term, - < w— Pv > , which appears in the governing equa- 

dy 

tion for the turbulent kinetic energy (Hinze, 1975). This term is some- 
times neglected, partly because it cannot be measured and partly because 
pressure tends to be poorly correlated with velocities, except near the 

wall (Townsend, 1956, and Tennekes and Lumley, 1972). Fig. 4.31 shows 

3 — 

the profile of the resolvable pressure work term, - < g^- Pv > . It can 
be seen that in the regions away from the wall (y > -.8), - < Pv > 

is much smaller than its corresponding values in the vicinity of the wall. 
In addition, the general shape of - < pv > is in accordance with the 
estimates of Laufer (1954) and Townsend (1956). These estimates were ob- 
tained from the turbulent energy balance in a pipe flow (see Chapter I). 
The average resolvable pressure velocity-gradient correlations 


3w 


< P t— > are shown 

az 


(pressure-strain terms), < P > , < P — > , and 

in Fig. 4.32. These terms govern the exchange of energy between the 

three components of resolvable turbulent kinetic energy. Note that since 

the sum of the above pressure velocity-gradient correlations is zero, 

these terms only transfer energy from one component to another, without 

changing- the total energy. Moreover, the negative sign for < P(3u./3x. > 

__ 2 l/2^ 

(no summation) indicates transfer of energy from < (u^ - < u^ >) > 

to other components (loss), whereas^a positive sign denotes energy gain. 

The profiles of < P > and < P g^- > show that throughout the channel 

the averaged streamwise component of resolvable turbulence intensity 

transfers energy to the other components, while the spanwise component 

receives energy. It is interesting to note that in the vicinity of the 

wall there is a large transfer of energy from the vertical component of 

turbulence intensity to the spanwise component. This is consistent with 

2 1/2 

the deficiency of the resolvable portion of the < v > profile in 
the region close to the wall shown in Fig. 4.28. 

In order to gain better insight into the flow of energy caused by 
the fluctuating pressure gradients, one might consider the governing 
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equations for each component of the resolvable turbulence energy. In 


these equations, the only terms where pressure appears explicitly are: 

— gp gp 

v r\ ^ . and ** ^ w — ^ for x_ v . v. compo - 


. 3P v. 


3y ' • ana - ^ w > tor x, y, 
nents of turbulence energy, respectively. Note that 


. - 3P . 

" < u 3^ > = 


< ? H> 


and 


. - 3P _ 
<w s?* 


<? u > 


but 


. - 3P . , 

<v & > * 


< ?| 2 > 

3y 


The average resolvable velocity pressure-gradient correlations are 
shown in Fig. 4.33. Examination of the < v > profile reveals that, 
aside from some energy loss in the region -.95 < y < -.83, the vertical 
component of the resolvable turbulent energy receives energy via - < ~v -p- > . 
Thus, - < v -■P’ > 


is primarily the source of energy for 


< v 2 > 1/2 . 


c/y 
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Chapter V 

CONCLUSIONS AND RECOMMENDATIONS 


In this work, we have numerically integrated the three-dimensional, 
time-dependent primitive equations of motion for the case of turbulent 
channel flow. To accomplish this task, a new, partially implicit algo- 
rithm and a new subgrid-scale model for' the inner region of the boundary 
layer were developed. An important feature of this partial implicit 
scheme is that the equation of continuity is solved directly. This, in 
turn, allows one to abandon the use of the Poisson equation for pressure. 
In addition, the stringent requirement on the time step caused by the 
numerical stability criterion for the diffusion equation is largely eased. 

The present computation has shown that many of the important fea- 
tures of wall-bounded turbulent flows can be reproduced using the Large 
Eddy Simulation approach. The overall agreement of the computed mean 
velocity and turbulence statistics with experimental data is satisfactory. 

In the present formulation of the subgrid scale model, the specifi- 
cation of the SGS length scale is not based on a well-def inted foundation. 
There are several choices available for this quantity which warrant sys- 
tematic study in this area; It would be desirable, for example, to in- 
corporate a Reynolds number dependence in the function defining the SGS 
length scale. This function, in turn, should allow for the vanishing of 
the subgrid scale model in a laminar flow. The profiles of total turbu- 
lent intensities indicate that, with the present grid resolution, a sub- 
grid scale model which allows anisotropy of SGS energy components is 
desirable. This modification of the subgrid scale model may not be nec- 
essary, if better grid resolution could be utilized. Nevertheless, the 
performance of the subgrid scale model used here is encouraging. 

In the light of our discussions about the grid resolution, a simula- 
tion with 32 x 65 x 128 mesh in x, y, and z direction, respectively, 
is strongly recommended. We believe that such a calculation will consid- 
erably improve the results obtained here and will provide the means for 
an objective evaluation as well as improvement of the subgrid scale model. 
It should be noted that this computation can presently be performed on 
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the ILLIAC IV computer. In addition, the use of a computer graphic sys- 
tem in conjunction with this simulation is highly desirable. This would 
provide the means for an efficient and a relatively convenient study of 
the detailed structures in the flow. Such a study, in turn, can consid- 
erably increase our knowledge of the structure and the mechanics of tur- 
bulent boundary layers . 

Based on the experience gained in our initial numerical experiment 
(Section 4.3), the following recommendations are made for the numerical 
simulation of laminar- turbulent flow transition: 

• Using an eddy viscosity model, the numerical simulation of transition 
from laminar to fully turbulent flow may be possible, provided that 
finite amplitude disturbances are added to the laminar flow. 

• However, if one wants to study the time evolution of small distur- 
bances, the eddy viscosity model should be used only after breakdown. 
Prior to breakdown, the use of any subgrid scale model may not be 
necessary. 

In extending the method to other flows, an important numerical prob- 
lem which must be resolved is the handling of inflow-outflow boundary 
conditions. In addition, an efficient numerical method should be devised 
which can be used in calculating flows that are inhomogeneous in more 
than one direction. Fully developed turbulent flow in a straight duct 
with a rectangular cross section is an example of such a flow. In simu- 
lating this flow, one can use periodic boundary conditions in the stream- 
wise direction. 

An important problem to study would be the case of turbulent flow 
over a smooth, flat plate. This flow is homogeneous only in one direc- 
tion. Moreover, its numerical simulation involves the handling of inflow- 
outflow boundary conditions. In addition, a suitable coordinate trans- 
formation should be used to map the infinite physical domain to a finite 
computational box. It is believed that the numerical simulation of this 
flow is an essential step towards the utilization of the Large Eddy Simu- 
lation approach in problems of engineering interest. 

It will be some time before the Large Eddy Simulation technique can 
be used in calculating flows of practical interest. However, in the in- 
terim, much information on the structure of turbulence can be obtained 
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by applying the method to simple but basic flows. The information, in 
turn, can be used in developing turbulence models in a simpler method 
for complex flows.' A knowledge of the pressure-velocity gradient corre- 
lations., for example, is of considerable value to the turbulence model- 
ers. As was shown in this study, using the Large Eddy Simulation ap- 
proach one can compute their large-scale components. Moreover, with the 
Large Eddy Simulation technique one can simultaneously obtain detailed 
quantitative information about the large-scale structures of the flow at 
thousands of spatial locations (grid points) throughout the flow field. 
This information cannot be gained from laboratory measurements . On the 
other hand, in the laboratory, one is capable of obtaining a long time 
history of the flow at relatively few spatial locations with minor ex- 
pense. With the present computers, this latter information about the 
flow can be gained only at high cost. Thus, at present, combined efforts 
of measurements and Large Eddy Simulation of turbulence seems to be an 
attractive approach to a better understanding of turbulent flows. 

The Large Eddy Simulation of turbulence is just beginning to emerge 
from its infancy, but it has already demonstrated a great potential in 
supplementing laboratory measurements. 
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Fig. 4.2. The resolvable portion of turbulence stress in the 
lower half of the channel at t = 0.45. 
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Fig. 4 
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.4. The resolvable portion of turbulence stress in the 
lower half of the channel at t = 0.85. 
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Fig. 4.8. Planar average of the resolvable portion of the stream- 
wise turbulence intensity in the vicinity of the lower 
wall at t = 1.05. 
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Fig. 4.9. Planar average of the resolvable portion of the stream- 
wise turbulence intensity in the vicinity of the lower 
wall at t = 1.425. 
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Fig, 4,10. Planar average of the resolvable portion of the stream- 
wise turbulence intensity in the vicinity of the lower 
wall at t = 2.025. 


62 











Fig. 4.12. Planar average of the resolvable portion of the span- 
wise turbulence intensity in the vicinity of the lower 
wall at 

A) t = 1.425 

B) t = 2.025 
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o Mean Velocity 

Instantaneous velocity profile at 

t = 1.625, x = 0 , z = I3h 3 

Instantaneous velocity profile at 

t = 1.825, x = 0 , z = 1 3h 3 
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Instantaneous streamwise velocity profiles obtained at one (x,z) 
location and at two different times. The instantaneous velocity 
profiles from the measurements of Grass (1971) are included in the 
left-hand side of the figure. 
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o Mean Velocity 

— Instantaneous velocity profile at 

t = 1.275, x = lOhi , z = I0h 3 

— Instantaneous velocity profile at 

t = l.05, x = IOfy, z=IOh 3 
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Fig. 4.15. Instantaneous streamwise velocity profiles obtained at one (x,z) 
location and at two different times. 




V 


Fig. 4.16. Instantaneous vertical velocity profiles obtained 
at the same location and times as in Fig. 4.15. 



Fig. 4.17. Instantaneous streamwise velocity profiles obtained at the same time 
and at two different (x,z) locations. 
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Fig. 4.22. Instantaneous spanwise variation of v (upper figure) and w (lower 
figure) at t = 1.05, y + = 3.85, x = 4h . 
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Fig. 4.23. Instantaneous spanwise variations of w at t = 1.05, x = 4h^, y = -.807 
(upper figure) and at t = 1.425, x = 4^, y = -.304 (lower figure). 





Fig. 4.25. Resolvable portion of turbulence stress averaged rn 
time. 
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Fig. 4.26a. Total turbulence stress (resolvable portion + SGS 
contribution), averaged m time. 
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4.26b. The resolvable and total turbulence stress in the 
vicinity of the walls. 

A) Near the lower wall 

B) Near the upper wall 
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Fig. 4.27. Mean velocity profile. The experimental data of Laufer, Comte-Bellot, and Hussain 
and Reynolds are included. 






Fig, 4.28, 


Time-averaged streamwise turbulence intensity in the vicinity of the wall (A) and 
away from the wall (B). 





Fig. 4.29. Time-averaged spanwise turbulence i 
away from the wall (B). 









Fig. 4.30. Time-averaged vertical component of turbulence intensity m the vicinity of the 
wall (A) and away from the wall (B) . 





Fig. 4.31. Pressure work term in the vicinity of the wall (A) 
and away from the wall (B). 
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Fig. 4.33. Velocity pressure-gradient correlations in the vicinity 
of the wall (A) and away from the wall (B) . 
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Appendix A 


FILTERING WITH NON-UNIFORM FILTER WIDTH 


In this appendix, we briefly discuss non-uniform width filtering 
and demonstrate its mathematical disadvantages. The use of such fil- 
ters (non-uniform width) is desirable in the directions in which the 
flow is inhomogeneous (see Section 2,1). For demonstration, we consider 
only simple box averaging as the filtering operation. 

Let 


f (x) = 


1 / *x+A,(x) 

(A + (x) + A_(x))/ x _ a (x) 


(A.1) 


where A + and A_ are the distances from x to its adjacent grid points. 
They will be treated as continuous variables. Note that here we consider 
only a one-dimensional case. Differentiating f yields: 
d 


9f 

9x 


d X < A + + A -> 


/ 


x+A, 


(A + (x) + A_(x)) J x-A 


ftt) d£ 


(A + + AJ 


/ 

dA \ / 

dA \ 

f(x + A + )^: 

+ dx ) - f <* - A ->( ] 

L "^f) 


or 


df 

dx 


(A + + A_) _ f (x + A + ) - f (x - A_) 

f + 


(A + A ) 

-r — 


(A + + AJ 


(A. 2) 


<A + + A _) 


dA. 


dA 


f(x + A.) — h f(x - A ) — — 

+ dx - dx 


Using the definition (A.l), we have: 


3f 

9x (A , 


^ /*x+A 

^La 


l + 9f 
3? 




A + + A 


f (x + A ) - f (x 

+ 


-A.)] 


(A. 3) 


Substitution of (A. 3) in (A. 2) yields: 
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3f 

3x 


St , jk (A + + A -> T 

8x (Aj + A ) 

T* — 


(A + + A_) 


dA 

f (x + A ) + f (x 

+ dx 


dA 1 

- o ttJ 


-CA.4-)- 


Thus., it -is clear that, in' general. 


3f , 3f 
3x f 3x ; 

The The above inequality and the presence of unfiltered quantities in 
(A. 4) renders the use of explicit nonuniform width-filtering extremely 
difficult. 
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Appendix B 


THE NUMERICAL DIFFICULTY WITH EXPLICIT TIME ADVANCING 
OF EQUATIONS OF MOTION 


In this appendix, we formally demonstrate the numerical difficulty 

associated with the fully explicit numerical integration of the Navier- 

* 

Stokes equations (see Section 3.3). Chebyshev polynomials and Fourier 
series are used to represent the flow variables in the vertical and 
horizontal difections, respectively. Consider the governing equations: 


3F , u 

u . = — -5 h H. 

x ox. x 

X 


(B.I) 


where H.. contains the transport and diffusion terms and a • over a 
variable denotes time derivative. Let 


P = 


i(k., x^k^x^) 
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n-=o k 7 

1 3 

th 


where T n (x 2 ) is the n -order Chebychev polynomial of the first kind 
and the double primes indicate that the first and last terms in the ser- 
ies are to be taken with factor 1/2. Eqns. (B.I), (B.2), (B.4), and 
(B.5) yield 


“x - S 

7. ? TT 


ik i a m (k i’ k 3^ 


-a ~(ki >k’ ) |+ C im VT m ( X2 ) e 


i(k iV k 3 X 3 ) 


m=o k l k 3 


m l 3 

r lk'a m (k ilk 'V 


(B.6) 


Other choices are possible. 
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From (B.3) we readily obtain: 


b in (k l’ k 3 ) “ N-N.N 2-f 2^ lu 

12 3 5=1 x 3 


H Z) “i (x r cos 0 j-’- x 3 ) 


• cos n0 . e 
3 


-i(k 3 _ Xl +k 3 x 3 ) 


(B. 7) 


where x n = cos 0. (see Eqn. (3.12)) and 0. = irj/N 0 , j = 0,1,2, ... ,N 0 . 

J 3 2 

Note that here we have enforced the no— slip boundary conditions, i.e. , 


u. (x. ,0. ,x,) 

J j=0,N 2 

Substituting (B.6) into (B. 7) , we get: 


= 0 


h 2 -i B 2 ( / llt l a m^ k i’ k p 

hJW - ssIEE -^ k i- k 3> 

1 2 3 j-1 Xl *3 m-o k 2 k 3 M 

. ' -lk 3 a m^ k l* k 3 

I _m« Mi 

+ C im | C0S C0S e 


The use of orthogonality of the expansion functions yields: 

/ -lk l a n ^1 ’ k 3 H ( / _:Lk l a m (k l ’ k 3 ^ ’ 


b in (k l’ k 3> = !- a A< k i> k 3 > 

\-ik 3 a n (ki, k 3 ) 


N 2 (/ 1 

' k £ jr^h’V ) 

-fk^a (V . k 


- ik 3 a m (k l’ k 3 ); 


l) nim + 1 


(B.9) 


The last term in (B.9), which is the result of enforcing the no-slip 
boundary conditions, is the source of trouble. To make this clear, con- 
sider the above equation for i = 1: 

b ln = - ik 1 a n (k 1 ,k 3 ) + C + a^,^) + (-l) n 3( k i> k 3 ) (B.10) 


a(k 1 ,k 3 ) “ Z ( " ik l a m^ k l ,k 3^ + C lm^ 
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N_ 

B( k i, k 3 ) = - ^ X) ^~ lk l a iEi^ k l’ k 3^ + 


Multiplying (B.10) by T n (cos 0..) and summing over all n yields: 


N2 


* k 3 ) — i k -^F ( > x 2 , k 3 ) "b H^( k -|^»x 2 , k 3 ) 

j 4 J 

+ a( k l5 k 3 ) cos n0 + 3(1^, k )^^ (-l) n cos n0. 

n=o 3 n=o 3 

where ^ over a variable denotes two-dimensional Fourier transform of 
that variable. But , 0, 


No 


sin H,6 S cot -f - 0 , j if 0 




n=o 


2 j 


~ 1 , 


(B .11) 


j - 0 


and 


0. 


No 


1 ' 

^ sin N 2 0 j tan = 0 , j ^ N 2 


y ' (-l) n cos n0. = 

n = o 1 


(B .12) 


n=o 


' N 2 ” 1 


j = No 


Note 


0 . 

J 


= Hi 


N„ 


j = 0,1,2,..., IL 


Hence, it has been shown that, unless 


01 ( 1 ^, Ic^) = 0 and BCk^,^) - 0 


(B.13) 


the two-dimensional discrete Fourier transform of u^ is discontinuous 
at the walls. It should be noted that (B.13) is equivalent to 


9P 

3x 


x 2~±l 


*1 


X2=±l 


which is the strearawise momentum equation evaluated at the walls (see 
Eqn. (3.15)). Similarly, the two-dimensional discrete Fourier transform 
of u^ or u^ is discontinuous unless 
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3P 

3y 


= H„ 


x 2 =±l 


I x 2 =±1 


or 


9P 

3z 


x 2 =±1 


H. 


x 2 =±1 


respectively. Therefore, if Neumann boundary condition is used for the 
Poisson equation, the Fourier transforms of u^ and will have dis- 

continuity at the boundaries. On the other hand, if Dirichlet boundary 
condition is used, the Fourier transform of u 2 will be discontinuous 
at the walls. In practice, the presence of discontinuity in the depen- 
dent variables results in non-convergent expansions which render a mea- 
ningless computation. A remedy for this problem is presented in Section 
3 . 4 . 
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Appendix C 


LISTING OF THE COMPUTER PROGRAM 
FOR THE CALCULATION OF TURBULENT CHANNEL FLOW 
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Cxxxxxxxxxxxxxxxxxxxx COMDECKS XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
CXXXXXXXXX BY CALLING A COMDECK THE DIMENSION OR COMMON STATEMENT * 

Cxxxxxxxxx FOLLOWING THE COMDECK WILL BE PLACED IN THE CALLING ROUTINEX 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

XCOMDECK Cl 

DIMENSION XB(16),YB(16.) 

XCOMDECK C2 

COMMON/AVEDY/RMIUC65) 

XCOMDECK C3 

C0MM0N/WV/WAVEX(16) ,WAVEYC16) »WAVEXS(16) ,WAVEYS(16) 

XCOMDECK C4 

DIMENSION BETA1(65),BETA2C65) 

XCOMDECK C5 

DIMENSION RHSV(4,61),AMB(4,4,61),AB(4,4,61),APB(4,4,61),AAUXC4,4, 
161),AMAUX(4,4,61),APAUX(4,4,61) 

XCOMDECK C6 

COMMON/S EC0ND/AP2( 65) , BP2C65) * CP2C65) 

XCOMDECK C7 

COMMON/L AGRNG/AP (65) , BP(65)»CP(65)»APR(65)»BPR(65)>CPR(65)»DPR(65) 
1 , EPRC65) 

XCOMDECK C8 

DIMENSION Z1(16,16),ZM1(16,16),D2<62) 

XCOMDECK C9 

DIMENSION BC1R(16, 16 ) » BC1I ( 16 » 16) »BCM1R(16 , 16 ) ,BCM1I ( 16 » 16 ) 
XCOMDECK CIO 

COMMON/DAT21/XR (16), XI (16) 

XCOMDECK Cll 

DIMENSION HR(16,16,65) 

LEVEL 2, HR 
XCOMDECK A1 

C0MM0N/DATA7/FR(16,16),FIC16,16) 

XCOMDECK A2 

COMMON DUDX(16,16,65) 

XCOMDECK A3 

CGMM0N/LCM4/DIVCQ6, 16 > 65) 

LEVEL 2, DIVC 
XCOMDECK A4 

COMMON/L ARGE2/PC 16, 16 ,65) 

LEVEL 2,P 
XCOMDECK A5 

COMMON/LARGEl/G( 16, 16,65) 

LEVEL 2, G 
XCOMDECK A6 

C0MM0N/LCM2/U(16,16,65),VC16,16,65),W(16,16,65) 

LEVEL 2, U , V, W 
XCOMDECK A 7 

COMMON/L CM1 /HI (16,16,65) ,H2(16 ,16,65),H3(16,16,65) 

LEVEL 2,H1,H2,H3 
XCOMDECK A8 

COMMON/LCM3/RU( 16 , 16,65),RV(16,16,65),RW(16,16,65) 

LEVEL 2, RU, RV, RW 
XCOMDECK A 9 

COMMON/STR/ZETA(65),Z(65)>RL(65),RM(65),E2,F2,EN,FN,R2,RN,A(65), 
1C(65),D(65), RR2 , RRN 
XCOMDECK A10 

DIMENSION G(16 ,16,65) 

LEVEL 2, G 
XCOMDECK All 

DIMENSION U1(16,16,65),U2U6,16,65),U3(16,16,65) 

LEVEL 2,U1,U2,U3 
XCOMDECK A12 

DIMENSION U(16,16,65),V(16,16,65),W(16,16,65) 

LEVEL 2,U,V,W 
XCOMDECK A13 

DIMENSION USUM(65),VSUM(65),WSUM(65) 

XCOMDECK B1 

C0MM0N/FLT/FILTX(16),FILTY(16) 

XCOMDECK B2 

COMMON/ EDDY/CV(63) 

XCOMDECK B3 
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COMMON/RECOVER/FACTOR( 65) 

XCOMDECK BA 

CQMMON/HGRI AV/U2S( 65) , V2S(65) »W2S(65),SSUM(65)» EDYVK65) 

XCOMDECK. B5 

COMMON/ AVEDY/MIUC65) 

XCOMDECK B6 

DIMENSION U2STC65) , V2ST (65 ) , W2ST(65) , UWT(65) 

XCOMDECK B7 

COMMON/SECOND/AP2(65),BP2(65),CP2C65) 

XCOMDECK B8 

COMMON/PENTA1/A1 (65) »B1(65)»C1(65) , Dl'(65) »E1C65),F1(65) 

*COMDECK B9 

DIMENSION U(16 , 16 , 65) 

LEVEL2,U 
KDECK MAIN 

PROGRAM MAIN ( INPUT, OUTPUT , TAPE8, TAPE9, TAPEIO ,TAPE1 1 ) 

CXXXXHTHIS SUBROUTINE MONITORS THE OVERALL SEQUENCE OF THE COMPUTATION 
C****X U,V,U ARE THE VELOCITIES IN STREAMWISE ,X, SPANWISE, Y, AND VERTICAL, 
2 DIRECTIONS. 

COMMON/LT A1/USUMC65) , UTSUMC65) , STSUMC65 ) ,U2SMT(65) , V2SMT (65) 

I , W2SMT (65),PVT(65),PUT(65), PUNST (65),PVNST(65), PWNST(65 ) , PUT (65) 
2»TCONT 

C0MM0N/LTA2/PDUT(65 ) , PDVT(65) ,PDWT(65) ,PDUNT(65), PDVNT (65) , PDWNT 
1(65) 

CQMMON/SGTT/5GST (65),ETED(65) ,U2STT(65) ,V2STT(65) ,W2STT(65) 
1,TSHGS, TSCNT 
COMMON/ COUNT/ I I CONT 
COMMON/SING/IMR, JMR, IMI , JMI 
COMMON/ADV/NTIME 

DIMENSION A3(61),B3(61),C3(61),D3(61),E3(61) 

COMMON/TINC/DT 

C0MM0N/PENTA2/XI,QI,GI,YI,QJ,GJ,XN,QIN,GIN,YN,QJN,GJN,Q2,Q3,RC1, 
1RC2, RPI , RP2, RP3, RPA 
*CALL Cl 
XCALL B5 

REAL MIU 

DIMENSION VAUX( A , 61 ) 

DIMENSION AX(3,3>61) , APX(3,3» 61) , AMX(3,3,61) , AXX(3,3,61) , 

1APXX(3 , 3,61), AMXX(3, 3, 61),VH(3,61) 

XCALL C3 
XCALL CA 
XCALL C5 
XCALL B7 
XCALL C7 

COMMON/ BC/CE1 , CE2, CE3 , CEA , CE5 , CE6 
COMMON/IDENTN/CODE 
XCALL A1 
XCALL A2 
XCALL A3 
XCALL AA 
XCALL A5 

COMMON/ CONST/CIO 0 ,C101,IJK,IJ,NHP1,HALF 
COMMON/ DATA 9/IMAX, JMAX, LMAX, NHALFX, NHAL FY 
C0MM0N/SCM2/LMAXP1 , D1 , D2 , D9 , DA , D5 , D6 
XCALL A6 
XCALL A7 
XCALL A8 
XCALL A 9 

C0MM0N/SCM3/DELTA1,DELTA2,RE,E 

INTEGER TIME, TEND 

TEND=200 

C0F=1.5 

DT=0 . 001 

NTIME=0 

CODE=2. 

CALL INITIAL 
CALL TRANS 
CC=1 ./( IMAXXJMAX) 

TP=0 . 5 
Cl=2. 0 
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C4=1.0 

LMAXM2=LMAX-2 

Lt1AXM3 = LMAX-3 

LHPl=LMAX/2+l 

ICONT=0 

LCONT=0 

LMAXM1=LMAX-1 

CALL INICON 

NTIME=1 

CALL INITIAL 

DO 300 TIME=1,TEND 

NTIME=TIME 

IC0NT=IC0NT+1 

IICONT=ICONT-20 

CALI CGURANT(DT,NTIME,TEND) 

CALL DIVG 
CALL RH5 

IF(NTIME.EQ.l) GO TO 360 
IF(IICONT.NE.O) GO TO 350 
ICQNT=0 
360 CONTINUE 
CALL STAT 
350 CONTINUE 

Cxxxxx DEFINE THE WAVE NUMBER INDEPENDENT ELEMENTS OF THE BLOCK - 
CXXXXXTRI DIAGONAL MATRIX 
DO 600 K=2,LMAX 
BETA1CK)=-C1/CDT*CE+MIUCK))) 

' BETA2CK)=-C1/CDTX(E+2.KMIU(K)I) 

600 CONTINUE 

CXXXXX DEFINE THE ELEMENTS OF THE TRIDIAGONAL MATRIX FOR THE CASE K1=K2=0. 
DO 800 K=1,LMAXM3 
KP2=K+2 

B3CK)=BP2CKP2)+BETA1CKP2) 

A3CK)=AP2CKP2) 

800 C3CK)=CP2CKP2) 

T=(Z(3)-ZC2))X0.5 

AK=1./CTP*DTXBETA2C3)) 

CXXXXX AMB,AB,APB ARE THE ELEMENTS OF THE BLOCK TRIDIAGONAL MATRIX. 

DO 640 M=l,4 
DO 640 N=1 , A 
DO 6 AO K = 1 > LMAXM3 
ABCM,N,K)=0. 

AMB(M,N,K)=0. 

APBCM,N,K>=0. 

6 AO CONTINUE 

DO 6 45 K=1 , LMAXM3 
KP2=K+2 

AB( 1 » I » K) =BP2CKP2)+BETA1(KP2) 

ABC2,2,K)=AB(1,1,K) 

AB(3,3,K)-BP(KP2) 

ABCA,A,K)=BETA2CKP2)XBP(KP2)XDTXTP 

AB(A,3,K)=BP2CKP2)+BETA2(KP2) 

645 CONTINUE 

AB(A,A,1)=CE2XBETA2(3)XDTXTP 
AB C 4 , 4 , LMAXM3 ) =CE5XBETA2 C LMAXM1 )*DT*TP 
AB(A,3,I)=BP2C3)+BETA2(3)XCE1 
AB(A,3,LMAXM3)=BP2(LMAXM1I+BETA2(LMAXMI)XCE6 
DO 650 K=l» LMAXM3 
KP2=K+2 

APB(1,1,K)=CP2(KP2> 

APB(2,2,K>=APB(1,1,K> 

APB(3,3,K)=CP(KP2) 

APBC4,4,K)=CPCKP2)XBETA2CKP2)*DT*TP 

APB(4,3,K)=CP2(KP2) 

AMB(1,1,K)=AP2CKP2) 

AMB(2>2,K)=AP2CKP2) 

AMB(4,3,K)=AP2CKP2) 

AMBCA,A,K)=AP(KP2)XBETA2(KP2)XDTXTP 

AMB(3,3,K)=APCKP2) 

650 CONTINUE 
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ORIGINAL’ PAGE IS 
CF poor QUALITY 

AMB ( 4 , 4 , LMAXM3 ) =C E4*B FT A2 ( LMAXM1 ) XDT*TP 
APBC4,4,1)=CE3*BETA2C3)XDT*TP 

CXXXX* DEFINE THE ELEMENTS OF THE, K1 = 0 , BLOCK TRIDIAGONAL SYSTEM 
DO 750 M=l,3 
DO 750 N=l,3 
DO 750 K=1 , LMAXM3 
AXCM,N,K)=0. 

APXCM,N,K)=0. 

AMXCM,N,K)=0. 

750 CONTINUE 

DO 752 K = l» LMAXM3 
AXC2,2,K)=ABC4,3,K) 

APXC2,2,K)=APBC4,3,K) 

AMXC2,2,K)=AMBC4,3,K> 

AXC3,1,K)=ABC2,2,K) 

APXC3,1,K)=APBC2,2,K) 

AMXC3,1,K)=AMB(2,2,K> 

AX (2, 3, K)=ABC4,4,K) 

APXC2,3,K3=APBC4,4,K) 

AMXC2,3,K)=AMBC4,4,K) 

AX(1,2,K)=AB(3, 3»K) 

APX(1,2,K)=APBC3,3,K) 

AMXC1,2,K)=AMBC3,3,K> 

752 CONTINUE 

C*xxxx DEFINE THE RHS OF THE BLOCK TRIDIAGONAL SYSTEM 
CALL VISCOSCU) 

DO 610 K=3, LMAXMI 
DO 610 J=1,JMAX 
DO 610 1=1 , IMAX 

Ua,J,K)=BETAl(K)x(U(I,J,K)+DTX(COFXHl(I,J,K)-0.5XRUa,J,K)))- 
1DUDXCI, J,K)XC4 
610 CONTINUE 

CALL VISCOSCV) 

DO 615 K=3, LMAXMI 
DO 615 J=1 , JMAX 
DO 615 1=1, IMAX 

VC I, J,K)=BETA1CK)*CVCI, J,K)+DTX(C0FXH2(I, J,K)-0.5*RVCI,J,K)))- 
1DUDXCI, J > K)XC4 
615 CONTINUE 

CALL VISCOSCW) 

DO 620 K=3, LMAXMI 
DO 620 J = 1 > JMAX 
DO 620 1=1, IMAX 

WCI,J,K)=BETA2CK)XCWCI,J,K)+DTXCCOFXH3CI,J,K)-0.5XRW(I, J,K)>)- 
1DUDXCI,J,K)XC4 
620 CONTINUE 

CXXXXX FOURRIER TRANSFORM 
DO 625 K=3, LMAXMI 
CALL MOV LEV CUC1,1,K),FR(1,1)»IJ) 

CALL FFTXC1.0) 

CALL FFTYCl. 0»CO 

CALL MQVLEVCFRC1,1),UC1,1,K),IJ) 

CALL MOVLEVCFICl,n,RU(l,l,K),IJ) 

CALL M0VLEV(VC1,1,K),FRC1,1),IJ) 

CALL FFTXC1.0) 

CALL FFTYCl. 0, CO 

CALL M0VLEVCFRC1,1),VC1,1,K),IJ) 

CALL MOVLEVCFI Cl , 1),RVC1,1,K),IJ) 

CALL M0VLEVCWC1,1,K),FRC1,1),IJ) 

CALL FFTXC1 . 0 ) 

CALL FFTYCl. 0, CO 

CALL MOVLEV C FRC 1 , 1) »W( 1 , 1 » K) » I J ) 

CALL M0VLEVCFIC1>1), RWC 1 , 1 , K) , I J ) 

625 CONTINUE 

Cxxxxxx THE REAL AND IMAGINARY PARTS OF THE FOURIER TRANSFORM OF THE RHS 

CXXXXXX OF THE BLOCK TRIDIAGONAL MATRIX IS COMPUTED . 

CfcXXXXX NOW DEFINE THE MATRIX ELEMENTS FOR EACH K1 AND K2 . 

NHP1X = I MAX/2+1 
NHP1Y=J MAX/2+1 
NHP2X=NHP1X+1 
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NHP2Y=NHP1Y+1 
DO 630 J=1,JMAX 
DO 630 I=1,IMAX 
WAV=WAVEXS(I)+WAVEYSCJ) 

IF(I.EQ.l.AND.J.EQ.l) GO TO 662 
IF(I.EQ.l) GO TO 410 
I-FCJ.EQ.l) GO TO 420 
IF( J .GE.NHP2Y) GO TO 430 
410 IF( J . LT . NHP2Y) GO TO 630 
GO TO 722 

420 IF( I . LT . NHP2X) GO TO 630 
GO TO 440 

430 IF(I.EQ.l.GR.I.EQ.NHPlX) GO TO 630 
440 CONTINUE 

DO 635 K=1 , LMAXM3 
KP2 = K+2 

CXXXXXX FIRST SOLVE FOR IMAGINARY PART OF U,V,AND REAL PART OF W AND P. 
RHSV(1,K)=RUCI,J,KP2) 

RHSV(2,K)=RV(I,J,KP2) 

RHSV(3,K)-0 . 

RHSVC4,K)=WCI, J,KP2) 

635 CONTINUE 

DO 647 K=1 t LMAXM3 
KP2=K+2 

AB ( 3 , 1 , K ) =-WA V EX Cl ) XIMAX 
ABC3,2,K)=-WAVEYCJ)XJMAX 
AB(1,4,K)=-AB(3,1»K)XBETA1CKP2)XDTXTP 

ABL2,4,tO=-AB(3,2,K)XBETAl(KP2)XDTXTP 

647 CONTINUE 

CXXXXX REARRANGING THE ROWS FOR CENTRAL DIFFERENCING 
DO 655 M=l,4 
DO 655 K=1 , LMAXM3 
AAUXC1,M,K)=ABC3,M,K) 

AAUX(4,M,K)=AB(1,M,K) 

AAUXC3,M,K)=ABC4,M>K) 

APAUX(1,M,K)=APB(3,M,K) 

APAUX(4,M,K)=APBC1,M,K) 

APAUXC3,M,K>=APB(4,M,K) 

AMAUXCl,M,K)=AMBC3,M,IO 

AMAUXC4,M,K)=AMBC1,M,K) 

AMAUXC3,M,K)=AMB(4,M,fO 

AAUX(2,M,K)=AB(2,M,K) 

AMAUX(2,M,K)=AMB(2,M,K) 

APAUX(2,M,K)=APBC2,M,K) 

655 CONTINUE 

DO 310 M=l,4 
DO 310 K=1 > LMAXM3 
310 VAUX(M> K)=RHSV (M> K) 

DO 315 K=1 j LMAXM3 
RHS V ( 1 » K) =V AUX C 3 » K ) 

RHSV(4,K)=VAUX(1,K) 

RHSVC3,K)=VAUX(4,K) 

315 CONTINUE 
IMR=I 
JMR=J 

CALL MTDAG C AMAUX, AAUX, APAUX, RHSV, 4, LMAXM3) 

DO 660 K=3 , LMAXM1 
KM2 = K-2 

RUa,J,K)=RHSV(l,KM2) 

RV(I,J,K)=RHSV(2,KM2) 

WCI,J,K)-RHSV(3,KM2) 

GCI,J,K)=RHSVC4,KM2) 

660 CONTINUE 

Cxxxx* COMPUTE THE REAL PART OF PRESSURE TRANSFORM AT THE WALL. 

G(I,J,2)=QI*G(I, J,3)+GIXG(I, J, 4)-C2.x(l .-CE 1 J/C AP( 3)xDT) )XW(I , J, 3 ) 
GCI,J,LMAX)=QINXGCI,J,LMAXM1)+GINXGCI,J,LMAXM2)-C2.*C1.-CE6)/ 
1(CPCLMAXM1)XDT))XW(I,J,LMAXM1) 

GO TO 630 
662 CONTINUE 

GCI, J,3)=TXW(I, J,3)*AK 
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GU,J,4) = C1.-TXBPC3))XGCI, J,3)/<T*CP(3) ) OF POOR QUALITY 

G(I,J,2)=0. 

G(I,J,1)=-CP(2)XGCI,J,3)/AP(2) 

DO 663 K=4,LMAX 
AK=1./CTPXDTXBETA2CK)) 

Ga,J,K+l)=(W(I, J,K)XAK-AP(K)XG(I,J,K-1)-BP(K)*GCI,J,K))/CP<K> 

663 CONTINUE 

DO 664 K=2,LMAX 

wa,j,K)=o. 

RU( I » J , K) =0 . 

RVCI,J,K)=Q. 

664 CONTINUE 
GO TO 630 

Cxxxxx SOLVE WHEN K1 = Q 
722 CONTINUE 

CX*XXX FIRST SOLVE FOR U, SIMPLE TRIDIAGONAL 
DO 724 K=1,LMAXM3 
724 D3(K)=RU(I, J,K+2) 

CALL TRIB(A3,B3,C3,E3,D3,tMAXM3) 

DO 726 K=3, LMAXM1 
726 RUCI,J,K)=D3(K-2) 

Cxxxxxx SOLVE FOR V,W,AND P, BLOCK TRIDIAGONAL 
DO 728 K=1 > LMAXM3 
KP2=K+2 
VH(1,K)=0. 

VH(2,K)=W(I,J,KP2) 

VHC3,K)=RVCI,J,KP2) 

728 CONTINUE 

DO 730 K=1,LMAXM3 
KP2=K+2 

AX(l,l,K)=-WAVEYCJ)xJMAX 
AX(3,3,K)=-AXa,l,K)XBETAl(KP2)XDTXTP 
730 CONTINUE 

DO 732 M=l,3 
DO 732 N=l,3 
DO 732 K=1 » LMAXM3 
AXX<M,N,K)=AXCM,N,K) 

APXXCM,N,K)=APXCM,N,K) 

AMXXCM,N,K)=AMX(M,N,K) 

732 CONTINUE 
IMR=I 
JMR=J 

CALL MTDAGCAMXX,AXX,APXX,.VH,3,LMAXM3> 

DO 734 K=3 » LMAXM1 
KM2=K-2 

RV(I»J»K)=VHC1> KM2) 

WCI,J,K)=VH(2,KM2) 

GCI,J,K)=VH(3,KM2> 

734 CONTINUE 

Cxxxxx COMPUTE THE REAL PART OF PRESSURE TRANSFORM AT THE WALL. 

GCI,J,2)=QIXG(I, J»3)+GIXG(I, J,4)-C2.xa.-CEli/(APC3)XDT))XWa, J,3) 
GCI,J,LMAX)=QINxg(I,U,LMAXM1)+GINXG(I,J,LMAXM2)-(2.X(1.-CE6)/ 

1 ( CP ( LMAXM1 ) *DT ) ) XWC I , J , LMAXMI ) 

630 CONTINUE 

DO 665 J =1 , JMAX 
DO 665 1=1 , IMAX 

IFCI. EQ.l.AND.J.EQ.l) GO TO 810 
WAV=WAVEXSCI)+WAVEYSCJ) 

IF(I.EQ.l) GO TO 510 

IFCJ. EQ.l) GO TO 520 

IF( J .GE.NHP2Y) GO TO 530 
510 I F( J . LI . NHP2Y) GO TO 665 
GO TO 736 

520 IFCI.LT.NHP2X) GO TO 665 
GO TO 540 

530 IFCI.EQ.l.OR.I.EQ. NHP1X) GO TO 665 
540 CONTINUE 

CXXXXXX NOW SOLVE FOR REAL PART OF U,V,AND IMAGINARY PART OF W AND P. 

DO 67 0 K=I j LMAXM3 
KP2=K+2 
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RHSV(1,K)=U(I,J,KP2) 

RHSV(2,K)=V(I,J,KP2) 

RHSV(3,K>=0. 

RHSV(4,K>=RWCI,J,KP2) 

670 CONTINUE 

DO 677 K=1 t LMAXM3 
KP2-K+2 

ABT3,1,K)=WAVEXCI)XIMAX 
AB(3,2,K>=WAVEYCJ)xJMAX 
AB(1,4,K)=-AB(3,1>K)XBETA1CKP2)XDTXTP 
ABC2,4,K)=-AB(3>2,K)XBETA1CKP2)XDTXTP 
677 CONTINUE 

exxxxx REARRANGING THE ROWS FOR CENTRAL DIFFERENCING 
DO 649 M=l,4 
DO 649 K=1 j LMAXM3 
AAUXC1,M,K)=AB(3,M,K) 

AAUXC4,M,K)=ABC1,M,K) 

AAUXC3,M,K)=AB(4,M,K) 

APAUXCl,M,K)=APB(3,M,IO 

APAUXC4,M,K)=APBC1,M,K) 

APAUX(3,M,K)=APB(4,M,K) 

AMAUXCl»Mf K)=AMB(3jM»K) 

AMAUXC4 ,M» K) = AMB Cl »M» K) 

AMAUXC3,M,K)=AMBC4,M,K) 

AAUXC2 > M,K)=AB(2,M,K) 

AMAUXC2,M,K)=AMBC2,M,K) 

APAUXC2>M>K)=APBC2>M,K) 

649 CONTINUE 

DO 320 M=l,4 
DO 320 K=1 > LMAXM3 
320 VAUXCM,K)=RHSVCM,IO 
DO 325 K = 1 > LMAXM3 
RHSVCl,K)=VAUXC3,fO 
RHSVC4,K)=VAUXC1>K) 

RHSVC3,fO=VAUX(4,K) 

325 CONTINUE 
IMI=I 
JMI“J 

CALL MTDAGCAMAUX,AAUX,APAUX,RHSV,4,LMAXM3) 

DO 690 K = 3> LMAXM1 
KM2=K-2 

Ua,J,K)=RHSVCl,KM2) 

VCI,J,K)=RHSVC2^KM2) 

RWCI,J,K)=RHSVC3,KM2) 

DUDXCI, J/ K)=RHSVC4,KM2 J 
690 CONTINUE 

CXXXXX COMPUTE THE IMAGINARY PART OF PRESSURE TRANSFORM AT THE WALL. 

DUDXCI,J,2)=QI*DUDX(I,J,3)+GI*DUDXCI,J,4)-C2.xa.-CE1)/CAPC3)XDT 

1>)XRWCI,J,3) 

DUDXCI,J,LMAX)=QINXDUDXCI,J,LMAXM1)+GINXDUDXCI,J,LMAXM2)-C2. 

IX ( 1 . -CE6 )/(CPC LMAXM1 )XDT) )XRW( I » J > LMAXM1) 

GO TO 665 

CXXXXX SIMPLE TRIDIAGONAL SOLUTION WHEN K1=0 AND K2=0. 

810 CONTINUE 

DO 820 K = 1 > LMAXM3 
820 D3(K)=UCIj J/K+2) 

CALL TRIBCA3,B3,C3,E3,D3,LMAXM3) 

DO 825 K=3,LMAXM1 
825 U(I,J,K)=D3(K-2) 

DO 830 K=1 » LMAXM3 
830 D3CK)=VCI,J,K+2) 

CALL TRIBCA3,B3,C3,E3,D3,LMAXM3) 

DO 835 K=3, LMAXM1 
835 VCI,J,K)=D3CK-2) 

GO TO 665 

CXXXXX SOLVE WHEN K1=C 
736 CONTINUE 

CXXXXX FIRST SOLVE FOR U, SIMPLE TRIDIAGONAL 
DO 738 K=1 > LMAXM3 
738 D3(K)=UCI>J>K+2) 


100 



CALL TRIBCA3,B3,C3,E3,D3,LMAXM3) 

DO 740 K=3, LMAXM1 
740 UCI,J,K)=D3CK-2) 

CXXXXXX SOLVE FOR V,W,AND P, BLOCK TRIDIAGONAL 
DO 742 K=1 , LMAXM3 
KP2=K+2 
VHC1,K)=0. 

VHC2,K)=RWCI» J,KP2) 

VH(3,K)=VCI,J,KP2) 

742 CONTINUE 

DO 744 K=1 , LMAXM3 
KP2=K+2 

AXC1, 1,K)=WAVEYCJ)XJMAX 
AXC3,3,K)=-AXC1,1,K)XBETA1CKP2)XDTXTP 
744 CONTINUE 

DO 746 M=l,3 

DO 746 N=1 , 3 

DO 746 K=1 > LMAXM3 

AXXCM,N,K)=AXCM,N,K) 

APXXCM,N,K)=APXCM,N,K) 

AMXXCM,N,K)=AMXCM,N,K) 

746 CONTINUE 
It1I=I 
JMI=J 

CALL MTDAGCAMXX,AXX,APXX,VH,3,LMAXM3) 

DO 748 K=3 , LMAXM1 
KM2=K~2 

VCI,J,K)=VHC1,KM2> 

RWCI,J,K)=VHC2,KM2) 

DUDXCI, J,K)=VH(3,KM2) 

748 CONTINUE 


ORIGINAL PAGE IS 
OE POOR QUALITY, 


DUDXCI, J, 2) =QIXDUDXCI,J,3)+GIXDUDXCI,J,4)-(2.XC1.-CE1)/CAPC3)XDT 
1 ) ) XRLK I , J , 3) 

DUDXCI, J,LMAX)=QINXDUDX(I,J,LMAXM1)+GINXDUDXC I , J>LMAXM2)-C2. 
lxa.-CE6)/CCPCLMAXMl)XDT))XRW(I, J,LMAXM1) 

665 CONTINUE 

DO 704 J=1 » UMAX 
DO 704 1=1 , IMAX 
WAV=WAVEXSCI)+WAVEYSCJ) 

IFCWAV.GT. 0.00001) GO TO 704 
DO 694 K=1 , LMAXP1 
RWCI,J,K)=0. 


694 DUDXCI, J,K)=0. 

704 CONTINUE 

Cxxxxx USE THE FACT THAT THE FLOW VARIABLES ARE REAL TO OBTAIN THE REM A I 
Cxxxxx -NING FOURIER COEFFICIENTS. 

DO 627 K=2 , LMAX 
DO 627 1=1, IMAX 
UC I , NHPIY, K) =0 . 

VC I , NHP1Y , K) =0 . 

WCI,NHP1Y,K)=0. 

RU C I , NHP1Y , K) =0 . 

RV C I , NHP1Y , K) =0 . 

RWCI,NHP1Y,K)=0. 

GC I , NHP1Y , K) =0 . 

DUDXCI, NHP1Y,K)=Q. 

627 CONTINUE 

DO 629 K=2, LMAX 
DO 629 J=1 , UMAX 
UCNHP1X, J,K)=0 . 

VCNHP1X, J,K)=0 . 

WCNHP1X, J,K)=0. 

RUCNHP1X, J,K)=0. 

RVCNHP1X,J,K)=0. 

RWCNHP1X, J,K)=0. 

GCNHP1X, J »K)=0 . 

DUDXCNHP1X, J , K)=0 . 

629 CONTINUE 

DO 550 K=2, LMAX 
DO 550 J=NHP2Y, JMAX 
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JJ=JMAX-J+2 

DO 550 I=NHP2X,IMAX 

II=IMAX-I+2 

UCII,JJ,K)=U(I, J,K) 

VCII,JJ,K)=VCI,J,K) 

well, JJ,K)=UCI,J,K) 

GC-I-I, JJ,K)=GCI,J,iO 

U(I,JJ,IO=U(II,J,K) 

VCI,JJ,IO=V(II,J,K) 

WCI,JJ,K)=WCII,J,lO 

GCI,JJ,K)=G(II,J,iO 

RUCII,JJ,K)=-RUCI,J,K) 

RVCII,JJ,K)=-RVCI,J,fO 

RWCII,JJ,K)=-RWCI, J,K> 

DUDXC II , J J >K)=-DUDXC I , J »K) 

RUCI, JJ,K)=-RUCII,J,K> 

RVCI, JJ,K)=-RVCIX,J,K) 

RWCI, JJ,K)=-RWCII, J,K> 

DUDXC I , J J ,K) =-DUDXC II »J » K) 

550 CONTINUE 

DO 560 K=2,LMAX 

DO 560 I=NHP2X,IMAX 

II=IMAX-I+2 

UCII,1»K)=UCI>1,K) 

VCII,1,K)=VCI,1,K) 

WCII,l,IO=WCI,l,IO 

GCII,1,K)=GCI,1,K) 

RUCII,1,K)=-RUCI,1,K) 

RVCII,l,K)=-RVCI,l,IO 

RWCII,l,K)=-RWCI,l,iO 

DUDXC II, 1,K)=-DUDXCI,1,K) 

560 CONTINUE 

DO 570 K=2,LMAX 

DO 570 J=NHP2X,JMAX 

JJ=JMAX~J+2 

UC1,JJ,K)=UC1,J,K) 

VCl,JJ,iO=VCl,J,IO 

WCl,JJ,IO-WCl,J,K) 

GC1,JJ,K5=GC1,J,K) 

RUC1,JJ,K)=-RUC1,J,K) 

RVCI, JJ,K)=-RVC1,J,K) 

RWC1, JJ,K)=-RWC1, J,K) 

DUDXC1, JJ,K) = “DUDX(1, J,K) 

570 CONTINUE 
CXXXXX INVERSE TRANSFORM 
DO 695 K=3 , LMAXM1 
CALL MOVLEVCUC 1 , 1 , K) > FRC1 , 13 > IJ7 
CALL MOVLEVCRUC1 , 1,K)> FI(I» 1) >IJ) 
CALL FFTXC-l - 0 ) 

CALL FFTYC-1 . 0 , CC) 

CALL MOVLEVCFRC 1 »l),UCl,l,fO,IJ) 
CALL MO V L EV C FI C 1 , 1 ) , RU C 1 , 1 > K ) > I J ) 
CALL MQVLEVCVC1,1,K),FRC1,1),IJ) 
CALL MOVL EVCRVC1,1,K),FIC1,1)»IJ) 
CALL FFTXC-l. 0) 

CALL FFTYC-1 . 0 , CC) 

CALL M0VLEVCFRC1,1),VC1,1,K),IJ) 
CALL M0VLEVCFI<1,1),RVC1,1>K),IJ> 
CALL M0VLEVCWC1,1,K),FRC1,1),IJ) 
CALL M0VLEVCRWC1,1,K),FIC1,1),IJ> 
CALL FFTXC-l. 0) 

CALL FFTYC-1. 0,CC) 

CALL M0VLEVCFRC1,1),WC1,1,K),IJ) 
CALL M0VLEVCFIC1,1),RWC1,1,K),IJ) 
695 CONTINUE 

DO 702 K = 1 , LMAXP1 
DO 703 J-l , JMAX 
DO 703 I=1,IMAX 
‘ FRCI, J)=GCI, J,K) 

FICI, J)=DUDXCI,J,K) 
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703 CONTINUE 

CALL FFTX(-l.O) 

CALL FFTYC-1 . 0 , CC) 

DO 705 J=1 , JMAX 
DO 705 1=1 , IMAX 
G(I,J,K)=FR(I,J) 

DUDXd, J,K) = FICI, J) 

705 CONTINUE 
702 CONTINUE 

Cxxxxxx STORE DATA CRU,RV,AND RW) FOR NEXT TIME STEP GRTGTjVAT re 

CALL PARTIALC1,G) WJSHjJiNAL rAGE ISs 

DO 710 K = 1 , LMAXP1 BE £OOR QUALITY/ 

DO 710 J=l» JMAX ' 

DO 710 1=1 , IMAX 

RUCI,J,K)=H1(I,J,K)+DUDXCI,J,K) 

710 CONTINUE 

CALL PARTIAL(2,G) 

DO 715 K=1 » LMAXP1 
DO 715 J=1,JMAX 
DO 715 1=1, IMAX 

RV(I,J,K)=H2CI,J,!Q+DUDXa,J,lO 
715 CONTINUE 

CALL PARTIALC3,G) 

DO 720 K=1 , LMAXP1 
DO 720 J=1 , JMAX 
DO 720 1=1, IMAX 

RUICI, J»K)=H3(I, J,K)+DUDX(I> J,K) 

720 CONTINUE 
CALL LTAVG 
LCONT=LCONT+l 
LLCONT=L CO NT-20 
IFCLLCONT.NE.G) GO TO 450 
LCONT=0 
CALL LTPR 
450 CONTINUE 
TP=0 . 5 
C4=l . 0 
Cl=2. 0 

200 FORMAT (1X,1P9E14.5) 

COF=l . 5 

CALL EXTERNC3, 1,R2,RR2) 

PRINT 400, TIME 
NHT=TEND/2 

IFCNTIME. EQ.NHT) CALL STAT 
400 FORMAT ( 3X, X TIME STEP=X,I3) 

IFCNTIME. NE. TEND) GO TO 300 
WRITEC9) U,V,W 

WRITE! 9) UTSUM, U2SMT, V2SMT , W2SMT , STSUM,PUT,PVT , PWT , PUNST, PVNST, 

1PUINST , SGST, ETED, U2STT, V2STT , W2STT , TCONT , TSHGS , TSCNT 
2, PDUT, PDVT, PDWT, PDUNT, PDVNT, PDWNT 
CALL STAT 
CALL LTPR 
STOP i 

300 CONTINUE 
STOP 
END 

XDECK PARTIAL 

SUBROUTINE PARTIAL CM, U> 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
C THIS SUBROUTINE COMPUTES THE PARTIAL DERIVATIVE OF U . M=1 CORRESPONDS 
C TO DERIVATIVE IN THE X-DIRECTION ,M=2 CORRESPONDS TO THE DERIVATIVE 
C IN THE Y-DIRECTION , AND M=3 CORRESPONDS TO THE DERIVATIVE IN THE Z-DIRECX 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

COMMON/IDENTN/CODE 

COMMON/DAT A9/IMAX, JMAX, LM AX, NHALFX, NHALFY 
COMMON/CONST/CIOO ,C1G1,IJK,IJ,NHP1,HALF 
XCALL A2 
XCALL A 9 
XCALL C7 
XCALL B9 


m 
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XCALL C3 
XCALL A1 

LMAXP1=LMAX+1 
DO 20 J=1 » JMAX 
DO 20 1=1 > IMAX 
DUDX(I,J,1)=0. 

DUDXC I , J , CMAXP1 ) =0 . 

20 CONTINUE 

IP CM.EQ.2) GO TO 30 
IF CM.EQ.3) GO TO 70 
Cxxxxxx DERIVATIVE IN THE X-DIRECTION 
DO 10 L=2» LMAX 
SIGN=1.0 

CALL M0VLEVCUC1,1,L),FR<1,1>,IJ) 

CALL FFTX(SIGN) 

DO 15 J=1 » JMAX 
DO 15 1=1, IMAX 
DUM=FICI,J) 

FICI, J)=WAVEX(I JXFRCI, J) 

FRCI,J)=-WAVEX(I)XDUM 
15 CONTINUE 
SIGN=-1.0 
CALL FFTX(SIGN) 

CALL MOVLEVCFR(l,l),DUDX(l,l,L),IJ) 

10 CONTINUE 
GO TO 300 
30 CONTINUE 

CXXXXXXDERIVATIVE IN THE Y-DIRECTION 
CC=1.0 

DO 35 L=2, LMAX 
SIGN=1.0 

CALL MOVLEVCUCl,l,L),FRa,l),IJ) 

DO 32 J=1 , JMAX 
DO 32 1=1, IMAX 
FICI, J)=0 . 0 
32 CONTINUE 

CALL FFTY(SIGN,CC) 

DO 40 J=l, JMAX 
DO 40 1=1, IMAX 
DUM=FICI,J) 

FI(I,J)=UAVEY(J1XFR<I,J) 

FR(I,J)=-WAVEYCJ)XDUM 
40 CONTINUE 
SIGN=-1.0 

CALL FFTYCSIGN, CC) 

CALL MOVLEVCFRU,l),DUDXtl,l,L),IJ> 

35 CONTINUE 
GO TO 300 
70 CONTINUE 

CXXXXXFIR5T DERIVATIVE IN THE Z-DIRECTION 
DO 82 J=1 , JMAX 
DO 82 1=1, IMAX 
DO 82 K=2, LMAX 
KP1=K+1 
KM1=K-1 

DUDXCI, J,K)=AP(K)XU(I , J,KM1)+CP(K)XU(I, J,KP1) 

82 CONTINUE 
90 CONTINUE 
300 CONTINUE 
RETURN 
END 

XDECK FFT 

IDENT FFT CA,B,N,ISN) 

ENTRY FFT 

X RADIX 2 COMPLEX FAST FOURIER TRANSFORM, COMPUTED IN PLACE. 

X SEE DON COMPUTING THE FAST FOURIER TRANSFORM, D R. SINGLETON, 
X COMM. ACM, V.10, N.10, PP. 647-654, OCT. 1967. 

X ARRAY A CONTAINS THE REAL COMPONENT OF THE DATA AND RESULT, 
x ARRAY B CONTAINS THE IMAGINARY COMPONENT. 

X N, THE NUMBER OF COMPLEX DATA VALUES, 


FFT2C 2 
FFT2C 3 
FFT2C 4 
FFT2C 5 
FFT2C 6 
FFT2C 7 
FFT2C 8 
FFT2C 9 
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'ORIGINAL PAGE I'.:. 
OS POOR QUALITY 


X 

MUST BE 

A POWER OF 2 AND 

GREATER THAN 1 


FFT2C 

10 

X 

THE 

SIGN OF ISN IS THE SIGN 

OF ' 

THE EXPONENTIAL IN THE TRANSFORM. 

FFT2C 

11 

X 

THE 

MAGNITUDE OF ISN IS 

THE 

INCREMENT SIZE FOR INDEXING 

FFT2C 

12 

X 

A 

AND B, 

AND IS ONE IN 

THE USUAL CASE. 


FFT2C 

13 

X 

DATA 

MAY ALTERNATIVELY BE STORED FORTRAN COMPLEX 


FFT2C 

14 

X 

IN 

A SINGLE ARRAY, IN 

WHICH CASE THE MAGNITUDE 


FFT2C 

15 

X 

OF 

ISN IS TWO AND ADDRESS 

B IS AC2), I.E. 


FFT2C 

16 

X 


CALL 

FFT2CA,AC2),N 

,2) 




FFT2C 

17 

X 

INSTEAD 

OF 





FFT2C 

18 

X 


CALL 

FFT2(A,B,N»1) 





FFT2C 

19 

X 

PROGRAM CONTAINS SINE TABLE 

FOR 

MAXIMUM N OF 32768 


FFT2C 

20 

X 

6400 

TIME 

FOR N=1024, 220 M 

.SEC 



FFT2C 

21 

X 

6400 

TIME 

FOR N=2XXM IS 

21. 

5XNXM MICRO-SEC. 


FFT2C 

22 

X 

6600 

TIME 

FOR N=1024, 44 

M. 

SEC. 



FFT2C 

23 

X 

6600 

TIME 

FOR N=2*XM IS 

4.3XNXM 

MICRO-SEC. 


FFT2C 

24 

X 

RMS 

ERROR 

FOR TRANSFORM- 

INVERSE 

IS LESS THAN 1.3E-13 


FFT2C 

25 

X 

FOR N=32768 OR SMALLER 





FFT2C 

26 

X 

FORTRAN 2. 

3 SUBROUTINE 





FFT2C 

27 

X 

BY R 

. C. SINGLETON, STANFORD RESEARCH INSTITUTE, NOV 

. 1968 

FFT2C 

28 


L100 

SXO 

B3 



NN 


FFT2C 

29 



SB4 

BO 



KK=0 


FFT2C 

30 



SB3 

B3-B7 



NN=NN-INC 


FFT2C 

31 



AXO 

1 



KSPAN=NN/2 


FFT2C 

32 



SB5 

BO 



K2=0 


FFT2C 

33 



SB6 

XO 





FFT2C 

34 



SXI 

B5 



K2”K2 


FFT2C 

35 



EQ 

B6 , B7 , FFT 



IFCKSPAN .EQ. INC) RETURN 

FFT2C 

36 


LUO 

SB4 

B3-B4 



KK=NN-KK 


FFT2C 

37 



SB5 

B3-B5 



K2=NN-K2 


FFT2C 

38 



SA2 

BI+B4 



EXCHANGE AtKK),ACK2) 

AND BCKK),BCK2) 

FFT2C 

39 



SA3 

B1+B5 





FFT2C 

40 



SA4 

B2+B4 





FFT2C 

41 



NX7 

X2 





FFT2C 

42 



SA5 

B2+B5 





FFT2C 

43 



NX6 

X3 





FFT2C 

44 



SA7 

A3 





FFT2C 

45 



SA6 

A2 





FFT2C 

46 



NX7 

X4 





FFT2C 

47 



NX6 

X5 





FFT2C 

48 



SA7 

A5 





FFT2C 

49 



SA6 

A4 



END OF EXCHANGE 


FFT2C 

50 



LT 

B6 , B4, LUO 



IFCKSPAN .LT. KK) GO 

TO LUO 

FFT2C 

51 

1120 

SB4 

B4+B7 



KK=KK+INC 


FFT2C 

52 



SB5 

B6+B5 



K2 = K5P AN+K2 


FFT2C 

53 



SA2 

B1+B4 



EXCHANGE ACKK),A(K2) 

AND B(KK),B(K2) 

FFT2C 

54 



SA3 

B1+B5 





FFT2C 

55 



SA4 

B2+B4 





FFT2C 

56 



NX7 

X2 





FFT2C 

57 



SA5 

B2+B5 





FFT2C 

58 



NX6 

X3 





FFT2C 

59 



SA7 

A3 





FFT2C 

60 



SA6 

A2 





FFT2C 

61 



NX7 

X4 





FFT2C 

62 



SXO 

B6 



K=KSPAN 


FFT2C 

63 



NX6 

X5 





FFT2C 

64 



SA7 

A5 





FFT2C 

65 



SA6 

A4 



END OF EXCHANGE 


FFT2C 

66 

L130 

AXO 

1 



K=K/2 


FFT2C 

67 



1X1 

X1-X0 



K2=K2-K 


FFT2C 

68 



PL 

XI , L130 



IFCK2 .GE. 0) GO TO 

L130 

FFT2C 

69 



LXO 

1 



K=K+K 


FFT2C 

70 



SB4 

B4+B7 



KK=KK+INC 


FFT2C 

71 



1X1 

X1+X0 



K2=K2+K 


FFT2C 

72 



SB5 

XI 



K2=K2 


FFT2C 

73 



GE 

B5,B4, LUO 



IFCK2 .GE. KK) GO TO 

LUO 

FFT2C 

74 



LT 

B4,B6,L120 



IFCKK .LT. KSPAN) GO 

TO L120 

FFT2C 

75 

FFT 







FFT2C 

76 



SB1 

XI 





INSR1 

1 



SA1 

Al+1 





INSR1 

2 



SB2 

XI 





INSR1 

3 
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L10 


L20 


130 


L40 

L50 


SAX 

A‘l+1 




INSR1 4 

SB3 

XI 




INSR1 5 

SAL 

Al+1 




INSR1 6 

SB4 

XI 




INSR1 7 

SA4 

B4 

ISN 



FFT2C 77 

MX2 

1 

MASK 



FFT2C 78 

SA5 

L60 




FFT2C 79 

SA3 

B3 

N 



FFT2C 80 

LX2 

57 




FFT2C 81 

PX7 

X3 




FFT2C 82 

BX6 

-X2*X5 




FFT2C 83 

PL 

X4» L10 

IFCISN .GE. 0) GO 

TO 

L10 

FFT2C 84 

BX6 

X2+X5 




FFT2C 85 

BX4 

-X4 

INC=-INC 



FFT2C 86 

LX3 

32 




FFT2C 87 

SA6 

A5 




FFT2C 88 

NXO 

B5,X3 




FFT2C 89 

PX2 

X4 




FFT2C 90 

SB7 

X4 




FFT2C 91 

DX7 

X2*X7 




FFT2C 92 

SA1 

B5+3 

SCM) 



FFT2C 93 

SB3 

X7 

NN=INC*N 



FFT2C 94 

SB6 

X7 

KSPAN=NN 



FFT2C 95 

EQ 

140 

GO TO L40 



FFT2C 96 

SA3 

CD 




FFT2C 97 

RX4 

X2*X1 

SDXCN 



FFT2C 98 

RX7 

X2XX0 

SDXSH 



FFT2C 99 

RX5 

X3XX0 

CDXSH 



FFT2C100 

RX6 

X3XX1 

CDXCN 



FFT2C1 0 1 

RX4 

X4-X5 




FFT2C102 

RX6 

X6+X7 




FFT2C103 

NX5 

X4 




FFT2C104 

RX7 

X1-X6 




FFT2C105 

RXO 

X0+X5 




FFT2C106 

NXX 

X7 




FFT2C1 07 

SB5 

B6+B4 

K2=KSPAN+KK 



F.FT2C1 08 

SA2 

B1+B4 

ACKK) 



FFT2C109 

SA3 

B1+B5 

A(K2) 



FFT2C110 

SA4 

B2 + B4 

BCKK) 



FFT2C111 

RX6 

X2+X3 




FFT2C112 

SA5 

B2+B5 

BCK2) 



FFT2C113 

RX2 

X2-X3 

RE 



FFT2C114 

SA6 

A2 

ACKK) 



FFT2C115 

RX7 

X4+X5 




FFT2C116 

RX3 

X1XX2 

CNXRE 



FFT2C117 

RX4 

X4-X5 

IM 



FFT2C118 

SA7 

A4 

BCKK) 



FFT2C119 

RX5 

X0XX4 

SNXIM 



FFT2C120 

RX2 

X0XX2 

SNXRE 



FFT2C121 

RX6 

X3-X5 




FFT2C122 

RX4 

X1XX4 

CNXIM 



FFT2C123 

SA6 

A3 

A(K2) 



FFT2C124 

RX7 

X2+X4 




FFT2C125 

SB4 

B6+B5 

KK=KSPAN+K2 



FFT2C126 

SA7 

A5 

BCK2) 



FFT2C127 

LT 

B4,B3,L30 

IFCKK .LT. NN) GO 

TO 

L30 

FFT2C128 

SB5 

B4-B3 

K2=KK-HN 



FFT2C129 

BX1 

-XI 

CN=-CN 



FFT2C130 

SB4 

B6-B5 

KK=KSPAN-K2 



FFT2C131 

LT 

B5,B4,L30 

IFCK2 .LT. KK) GO 

TO 

L30 

FFT2C132 

SB4 

B4+B7 

KK=KK+INC 



FFT2C133 

SA2 

SD 




FFT2C134 

LT 

B4,B5,L20 

IFCKK .LT. K2) GO 

TO 

L20 

FFT2C135 

SB4 

BO 

KK=Q 



FFT2C136 

SX5 

B6 




FFT2C137 

AX5 

1 

KSPAN=KSPAN/2 



FFT2C138 

SB6 

X5 




FFT2C139 

SB5 

B6+B4 

K2=KSPAN+KK 



FFT2C140 

SA2 

B1+B4 

ACKK) 



FFT2C1 41 

SA3 

B1+B5 

ACK2) 



FFT2C142 
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j^AGE S 

QE P-QOR 





SA4 

B2+B4 

B(KK) 

FFT2C143 


RX6 

X2+X3 


FFT2C144 


S AS 

B2+B5 

B(K2) 

FFT2CI45 


RX7 

X2-X3 


FFT2C146 


SA6 

A2 

ACKK) 

FFT2C147 


SA7 

A3 

A(K2) 

FFT2C148 


RX6 

X4+X5 


FFT2C149 


SB4 

B6+B5 

KK=KSPAN+K2 

FFT2C150 


RX7 

X4-X5 


FFT2C151 


SA6 

A4 

B(KK) 

FFT2C152 


SA7 

AS 

B(K2) 

FFT2C153 


LT 

B4,B3,L50 

IF( KK .LT. NN) GO TO L50 

FFT2C154 


EQ 

B6 , B7 , L10 0 

IFCKSPAN .EQ. INC) GO TO L100 

FFT2C155 


SA1 

A1 

SCM) 

FFT2C156 


SB4 

B7 

KK=INC 

FFT2C157 


RX6 

X1XX1 


FFT2C158 


$A1 

Al+1 

M=M+1, S(M) 

FFT2C159 


FX6 

X6+X6 


FFT2C160 


SA3 

ONE 


FFT2C161 


SA6 

CD 

CD=2XSCM)XX2 

FFT2C162 

L60 

BXO 

XI 

SN=SD 

FFT2C163 


RX6 

X3-X6 

CN=1 , 0-CD 

FFT2C164 


BX7 

xo 


FFT2C165 


NX1 

X6 


FFT2C166 


SA7 

SD 


FFT2C167 


EQ 

L30 

GO TO L30 

FFT2C168 

S 

DATA 

9.587 37 990 95977346 E~5 

FFT2C16 9 


DATA 

1 . 917 4759731 07 0331 E-4 

FFT2C17 0 


DATA 

3.8349518757139559E-4 

FFT2C171 


DATA 

7.6699031874270453E-4 

FFT2C172 


DATA 

1 .533 98 0186 2847 656 E~3 

FFT2C173 


DATA 

3. 067 956762965 9 763E-3 

FFT2C174 


DATA 

6. 135 88 464915447 54 E“3 

FFT2C175 


DATA 

1. 227153828571 9926 E-2 

FFT2C176 


DATA 

2.4541228522912288E-2 

FFT2C177 


DATA 

4. 9067 6743274180 14 E-2 

FFT2C178 


DATA 

9. 80 17 14032956 06 02E-2 

FFT2C17 9 


DATA 

1. 950 90 32201612827 E“1 

FFT2C180 


DATA 

3 . 826 8 34323 6 5 08 977 E“1 

FFT2C181 


DATA 

0.7071067811865475 

FFT2C182 

ONE 

DATA 

1.0 


FFT2C183 

CD 




FFT2C184 

SD 




FFT2C185 


END 



FFT2C186 

XDECK 

FFTX 





SUBROUTINE 

FFTXCSIGN) 



CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

C FAST 

FOURIER TRANSFORM IN 

X-DIRECTION 

X 


CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
C0MM0N/DATA9/IMAX, JMAX, LMAX»NHALFX,NHALFY 
XCALl A1 
XCALL CIO 

ISN=-SIGN 

IF (SIGN .LT. 0.) GO TO 3 
DO 2 J=1,JMAX 
DO 1 1=1, IMAX 
FI(I,J)=0. 

1 CONTINUE 

2 CONTINUE 

3 CONTINUE 

DO 100 J = 1 , JMAX 
DO 110 1=1 , IMAX 
XR(I)=FR(I, J) 

XI(I)=FI(I, J) 

110 CONTINUE 

CALL FFTCXR,XI,IMAX,ISN) 

DO 120 1=1, IMAX 
FR( I , J ) =XRC I ) 

Fia,j)=xicn 

120 CONTINUE 
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100 CONTINUE 
RETURN 
END 

XDECK FFTY 

SUBROUTINE FFTYCSIGN, COEF3) 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
C FAST FOURIER TRANSFORM IN Y-DIRECtlON X 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
XCALL A1 
XCALL CIO 

COMMON/DATA 9/IMAX, JMAX, LMAX,NHALFX,NHALFY 
ISN=-SIGN 
C Y-TRANSFORM 

DO 100 1=1 , IMAX 
DO 110 J=1,JMAX 
XR( J ) =FR( I , J ) 

XICJ)=FI(I, J) 

1-10 CONTINUE 

CALL FFT(XR,XI, JMAX,ISN) 

IFCSIGN.LT. 0.) GO TO 200 
DO 120 J=1»JMAX 
FRCI,J)=XRCJ) 

FICI, J)=XICJ) 

120 CONTINUE 
GO TO 100 

200 DO 130 J=1,JMAX 

FRC I > J ) =XRC J JXC0EF3 
FICi;j)=XICJ)XCOEF3 
130 CONTINUE 
100 CONTINUE 
RETURN 
END 

XDECK INITIAL 

SUBROUTINE INITIAL 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
cx THIS SUBROUTINE COMPUTES THE VARIOUS NECESSARY ARRAYS AND CONSTANTS X 
CXFOR SGS, PARTIAL, POISON, AND FILTER SUBROUTINES X 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
XCALL A1 

COMMON/ADV/NTIME 
XCALL B1 

COMMON/ DATA9/IMAX, JMAX, LMAX»NHALFX,NHALFY 
C0MM0N/SCM2/LMAXP1 , D1,D2,D3,D4, D5, D6 
C0MM0N/SCM3/DELTA1, DELTA2, RE, E 
C0MMQN/SCM4/CI , CJ , CK, CJK, CIK, Cl J 
XCALL C3 

COMMON/CONST/Cl 00,C101 ,IJK, IJ,NHP1»HALF 

REAL NAVG 

C=0.A 

S=2./3. 

PAI=AC0S<-1.) 

Cxxxxx DELTA1 AND DELTA2 ARE THE MESH SIZES IN X AND Y DIRECTIONS 
DELTAl=PAI/8. 

DELTA2=PAI/12 . 

IMAX=16 

JMAX=16 

LMAX=6A 

IJ=IMAXXJMAX 

LMAXP1=LMAX+1 

IJK=IMAXXJMAXXLMAXP1 

CI=1./IMAX 

CJ=1./JMAX 

CK=1 ./LMAXP1 

CJK=1./CJMAXXLMAXP1) 

CIK=1./CIMAXXLMAXP1) 

CIJ=1./CIMAXXJMAX) 

RE=6A0 . 25 
E=1 ./RE 
NHALFX=IMAX/2 
NHALFY=JMAX/2 
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NHP1X=NHALFX+1 
NHP1Y=NHALFY+1 
C100=2.GxPAI/CIMAXXDELTAl) 

Cl 0=2. OXPAI/C JMAXXDELTA2) 

C101=C100/IMAX 
C11=C1 0/ JMAX 

Cxxxxx DEFINE WAVE NUMBERS. 

Cxxxx* NOTE THAT UAVEX AND WAVEY ARE SMALLER THAN THE ACTUAL WAVE NUMBERS 
CXXXXX BY FACTOR OF IMAX AND JMAX RESPECTIVELY. 

DO 100 1=1, IMAX 
MM=I/NHP1X 


M=MMX IMAX+1 

WAV EX ( I ) =C101X( I-M) 

WAVEXSCI)=CC1Q0X(I-M))XX2 
100 CONTINUE 

WAVEX(NHPIX) =0 . 


WAVEXSCNHP1X) =0 . 

DO 130 J=1 , JMAX 
MM=J/NHP1Y 
M=MMXJMAX+1 
WAVEYC J ) =C11X( J-M) 

WAVEYSC J)=(C10X( J-M))XX2 
130 CONTINUE 

WAVEYCNHP1Y) =0 . 

WAVEYS(NHP1Y)=0. 

1000 FORMATC1P8E15.7) 

NAVG=2 

IF(NTIME.EQ.O) NAVG=6 

NHP2X=NHP1X+1 

NHP2Y=NHP1Y+1 

CXXXXXCOMPUTE THE NORMALIZED FOURIER TRANSFORM OF THE FILTER FUNCTION IN X-DIREC 
DO 300 J=1 » JMAX 
DO 300 I=1,NHP1X 

FR(I,J3=EXPC-6.XFLOAT(I-l)XX2/CNAVGXX2)) 

300 CONTINUE 

DO 310 J=l, JMAX 
DO 310 I=NHP2X,IMAX 
II=IMAX-I+2 
FRCI,J)=FRCII,J) 

310 CONTINUE 

CXXXXX COMPUTE THE NORMALIZATION CONST, AREA. 

AREA=0. 

DO 320 1=1, IMAX 
AREA=AREA+FR( I ,1 ) 

320 CONTINUE 

DO 330 J=1 , JMAX 
DO 330 1=1, IMAX 
FRCI, J)=FR(I, JI/AREA 
FI(I,J)=0. 

330 CONTINUE 

CALL FFTXC1.0) 

DO 340 1=1, IMAX 
FILTXCI)=FRCI,1> 

340 CONTINUE 

CXXXXXCOMPUTE THE NORMALIZED FOURIER TRANSFORM OF THE FILTER FUNCTION IN Y-DIREC 
DO 400 J=1 , NHP1Y 
DO 400 1=1, IMAX 

FR( I , J ) =EXP C “6 . X FLOAT ( J-l JXX2/CNAVGXX2) ) 

400 CONTINUE 

DO 410 J=NHP2Y, JMAX 
DO 410 1=1, IMAX 
JJ=JM AX-J+2 
FRCI, J)=FRCI, JJ) 

410 CONTINUE 
AREA=0 . 

DO 420 J=1 , JMAX 
AREA=AREA+FRC1,J) 

420 CONTINUE 

DO 430 J=1 , JMAX 
DO 430 1=1, IMAX 
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FRCI, J)=FRCI,J)/AREA 
FICI, J>=0. 

430 CONTINUE 

CALL FFTYCI . 0,1.0) 

DO 440 J=1,JMAX 
FILTY(J)=FR(1,J) 

440 CONT-INUE 

FILTX(NHP1X)=0. 

FI LTY(NHPIY) =0 . 

PRINT 1000, (WAVEXCL),L=1,IMAX) 

PRINT 1000, CWAVEY(L),L=1,JMAX) 

PRINT 1000, (WAVEXS(L),L=1,IMAX> 

PRINT 1000, CWAVEYSCL),L=1,JMAX) 

RETURN 

END 

XDECK INICON 

SUBROUTINE INICON 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

CX THIS SUBROUTINE GENERATES THE INITIAL FIELD FOR THE COMPUTATION X 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

C0MM0N/DATA9/IMAX, UMAX, LMAX, NHALFX»NHALFY 
DIMENSION GC161),Y(161),F(65) 

COMMON/SCM3/DELTA1, DELTA2»RE»E 
XCALL A8 

C0MM0N/CQNST/C100, C101 ,IJK,IJ,NHP1,HALF 
XCALL A1 
XCALL A13 

COMMON/SCM4/CI , CJ, CK, CJK, CIK,CI J 
XCALL All 

EQUIVALENCE CUI , HI ) , (U2,H2) , CU3 ,H3) 

XCALL A2 
XCALL A6 
XCALL A7 
XCALL A 9 

PAI=ACOS(-l.) 

LMAXP1=LMAX+1 
LMAXM1=LMAX-1 
DO 210 J=1 , JMAX 
DO 210 1=1 , IMAX 
U1CI,J,2)=0. 

U2CI, J,2) = 0. 

U3CI,J,2)=0. 

U1CI,J,1>=0. 

U2CI, J,1)=0. 

U3CI, J, 1 ) =0 . 

U1CI,J,LMAX)=Q. 

U2(I, J,LMAX)=0 . 

U3CI, J , LMAX) =0 . 

U1 C I, J , LMAXP1 ) = 0 . 

U2CI, J,LMAXP1)=0 . 

U3(I,J,LMAXP1)=0. 

UCI, J,2)=Q. 

V(I,J,2)=0. 

WCI,J,2)=0. 

U(I,J,1)=Q. 

VCI,J,1)=0. 

WCI,J,1)=0. 

U(I, J , LMAX) =0 . 

V(I,J,LMAX)=0. 

W(I,J,LMAX)=0. 

UCI, J, LMAXP1) = 0 . 

' VCI,J,LMAXP1)=0. 

W(I,J,LMAXP1)=0. 

210 CONTINUE 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

cx THE VELOCITY FIELD FOR THE INITIATION OF THE PROGRAM IS OBTAINED x 

CX FROM THE DISK. THE ORIGINAL VELOCITY FIELD IS GENERATED FROM A X 

CX SEPARATE PROGRAM (SEE SECTION 4.2 IN THE TEXT). X 

CX UI,U2,U3 ARE THE COMPONENTS OF THE VELOCITY FIELD AT TIME STEP N x 

Cx RU,RV, AND RW ARE THE INFORMATION AT TIME STEP N-l, NECESSARY FOR X 
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■ORIGINAL PAQB i $ 

OEEOOR QIJAUW 

CX ADAMS BASHFORTH METHOD. X 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
READ(8) U1,U2,U3,RU,RV,RU 
DO 25 K=2,LMAX 
DO 25 J-l , JMAX 
DO 25 1=1 , IMAX 
U(I,J,K)=U1(I,J,K) 

VCI, J,K)=U2(I, J,fO 
W(I, J,K)=U3CI,J,K) 

25 CONTINUE 

CALL EXTERN(3,1,R2,RR2) 

CALL EXTERNC31,33,RN,RRN) 

PRINT 2000 

1000 FORMAT ( 1P8E15 . 7 ) 

2000 F0RMAT(1H1,X VELOCITY IN THE X-DIRECTION ACCROSS THE CHANNEL X) 
PRINT 1000,(U(10,10,K),K=1,LMAXP1) 

RETURN 

END 

XDECK CURL 

SUBROUTINE CURLCU,V,W) 

Cxxxxx THIS SUBROUTINE COMPUTES THE VQRTICITY FIELD 
C0MM0N/DATA9/IMAX, JMAX, LMAX, NHALFX, NHALFY 
COMMON/ CONST/C1 00,C101,IJK,IJ,NHP1,HALF 
XCALL All 

EQUIVALENCE (Ul, HI ) , (U2, H2) , (U3, H3) 

XCALL A12 
XCALL A7 
XCALL A2 

LMAXP1=LMAX+1 
CALL PARTIAL (2, W) 

CALL MOVLEVCDUDXCl,l,l),Uia,l,l),IJJO 
CALL PARTIALC3, V) 

DO 10 K=1 , LMAXP1 
DO 10 J=1 , JMAX 
DO 10 1=1, IMAX 

UKI, J,K)=U1(I,J,K)“DUDX(I, J,K) 

10 CONTINUE 

CALL PARTIAL(3, U) 

CALL MQVLEV(DUDXC1,1,1>,U2(1,1,1),IJK) 

CALL PARTIAL ( 1 , W) 

DO 15 K=1 , LMAXP1 
DO 15 J=1 , JMAX 
DO 15 1=1, IMAX 

U2(I, J,K)=U2CI,J,K)-DUDX(I,J,K) 

15 CONTINUE 

CALL PARTIAL(1,V) 

CALL MOVLEVC DUDX( 1 ,1,1)>U3(1,1,1),IJK) 

CALL PARTIAL(2,U) 

DO 20 K=1 , LMAXP1 
DO 20 J=1 , JMAX 
DO 20 1=1, IMAX 

U3CI,J,K)=U3CI, J,K)-DUDX(I,J,K) 

20 CONTINUE 
RETURN 
END 

XDECK RHS 

SUBROUTINE RHS 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX^XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
Cx THIS SUBROUTINE COMPUTES THE RIGHT HAND SIDE OF THE GOVERNING 
CX EQUATIONS, EXCLUDING THE PRESSURE. 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

C0MM0N/DATA9/IMAX, JMAX, LMAX, NHALFX, NHALFY 
COMMON/CONST/CIO 0 , CIO 1 , 1 JK , I J , NHP1 , HALF 
C0MM0N/SCM2/LMAXP1 , D1 , D2 , D3 , D4 , D5 , D6 
C0MM0N/SCM3/DELTA1 , DELTA2, RE, E 
XCALL A2 
XCALL A5 
XCALL A6 
XCALL A 7 

CALL SGS 
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X * 



CXXX*XMOMENTUM EQUATION IN THE X-DIRECTION 
CALL PARTIAL d , V) 

DO 10 K = 1 , LMAXPI 
DO 10 J = 1 > UMAX 
DO 10 1=1, IMAX 

Gd,J,K)=Vd,J,K)XDUDXd,J,IO 
10 CONTINUE 
C 

CALL PARTIALCl.U) 

DO 20 K=l, LMAXPI 
DO 20 J=l> JMAX 
DO 20 I=1,IMAX 

Gd,J,K)=Gd,J,K)+Wd,J,K)*DUDXd,J,K) 

20 CONTINUE 

CALL PARTIALC2,U) 

DO 30 K = l, LMAXPI 
DO 30 J = 1 » UMAX 
DO 30 I=1,IMAX 

Gd,J>K)=Gd,J,IO-Vd,J,K)XDUDXd,J,K) 

30 CONTINUE 

CALL PARTIALC3,U) 

DO AO K = l, LMAXPI 
DO AO J =1 , UMAX 
DO AO 1 = 1 , IMAX 

Gd,J,K)=Gd,J,K)-Wd, J,K)XDUDXd,J,K> 

AO CONTINUE 

CALL FILTER(G) 

DO A5 K=l, LMAXPI 
DO A5 J=l, JMAX 
DO A5 1=1, IMAX 

Hld,J,K)=Gd, J,K)+HHI,J,K)+1. 

A5 CONTINUE 

CXXXXX COMPUTE THE VISCOUS TERMS IN THE X-MOMENTUM EQUATION 
CALL PARTIAL ( 1 , U) 

CALL M0VLEV(DUDXC1,1,1),G(1,1,1),IJK) 

CALL PARTIAL C 1 , G) 

DO 50 K=l, LMAXPI 
DO 50 J=1 , JMAX 
DO 50 1=1, IMAX 

Hld,J,K)=Hld, J,K)+EXDUDX(I,J,K) 

50 CONTINUE 

CALL PARTIAL(2»U) 

CALL MOVLEVCDUDXCl,l,l),Gd,l,l),IJK) 

CALL PARTIAL62, G) 

DO 55 K=l, LMAXPI 
DO 55 J = 1 > JMAX 
DO 55 1=1, IMAX 

Hid, J,K)=H1(I,J,K)+EXDUDXCI,J,K) 

55 CONTINUE 

CXXXXXMOMENTUM EQUATION IN THE Y-DIRECTIQN 
CALL PARTIAL(2,U) 

DO 65 K=l, LMAXPI 
DO 65 J=l, JMAX 
DO 65 1=1, IMAX 

G(I,J,K)=Ud,J,K)XDUDX(I,J,K) 

65 CONTINUE 

CALL PARTIAL(2,W) 

DO 70 K=l, LMAXPI 
DO 70 J=1 , JMAX 
DO 70 1=1, IMAX 

Gd,J,K)=Gd,J,K)+Wd,J,K)XDUDXd,J,K) 

70 CONTINUE 

CALL PARTIALC3, V) 

DO 75 K=l, LMAXPI 
DO 75 J=i, JMAX 
DO 75 1=1, IMAX 

Gd, J,K)=GCI, J,K)-WCI,J,K)»DUDX(I, J,K) 

75 CONTINUE 

CALL PARTIAL C 1 , V ) 

DO 80 K=l, LMAXPI 
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DO 80 J=1,JMAX 
DO 80 1=1 , IMAX 

GCI,J,K)=GCI,J,K)-UCI,J,K)XDUDXU, J,K) 

80 CONTINUE 

CALL FILTER(G) 

DO 85 K=1,LMAXP1 

DO 85 J=1 , JMAX 

DO 85 1=1, IMAX 

H2(I, J,K)=H2CI,J,K)+GCI, J,K) 

85 CONTINUE 

CXXXXXCOMPUTE THE VISCOUS TERMS IN THE Y-MOMENTUM EQUATION 
CALL PARTIALC1,V) 

CALL MOVLEVCDUDX(l,l,l),GU,l,l),IJK) 

CALL PARTIALU,G) 

DO 90 K=1 » LMAXP1 
DO 90 J=1 , JMAX 
DO 90 1=1, IMAX 

H2CI, J,K)=H2(I,J,K)+EXDUDXCI,J,K) 

90 CONTINUE 

CALL PARTIAL(2,V) 

CALL MOVLEVCDUDXa,l,l),G(l,l,l),IJK) 

CALL PARTIALC2,G) 

DO 95 K=1 , LMAXP1 
DO 95 J=1,JMAX 
DO 95 1=1, IMAX 

H2CI, J,K)=H2(I, J,K)+EXDUDXCI,J,K) 

95 CONTINUE 

CxxxxxMOMENTUM EQUATION IN THE Z-DIRECTION 
CALL PARTIALC3, V) 

DO 105 K = 1 , LMAXP1 

DO 105 J=1,JMAX 

DO 105 1=1, IMAX 

G(I, J,K)=VCI, J,K)XDUDXCI, J,K) 

105 CONTINUE 

CALL PARTIALC3,U) 

DO 110 K=1 , LMAXP1 
DO 110 J=1,JMAX 
DO 110 1=1 , IMAX 

GCI,J,K)=G(I,J,K)+UCI, J,K)XDUDXCI,J,K) 

110 CONTINUE 

CALL PARTIAL (2, W) 

DO 115 K=1 , LMAXP1 
DO 115 J=1 > JMAX 
DO 115 1=1, IMAX 

GCI, J,K)=G(I» J,K)-V(I, J,K)XDUDX(I, J»K) 

115 CONTINUE 

CALL PARTI AL ( 1 , U) 

DO 120 K=1 , LMAXP1 
DO 120 J=1 , JMAX 
DO 120 1=1, IMAX 

GCI, J,K)=GCX,J,K)-UCI,J,K)XDUDXCI, J,K) 

120 CONTINUE 

CALL FILTERCG) 

DO 125 K=1 , LMAXP1 
DO 125 J=1 , JMAX 




DO 125 1=1, IMAX 
H3(I»J,K)=H3CI,J»K)+GCI,J,K) 

125 CONTINUE 

CXXXXXCOMPUTE THE VISCOUS TERMS IN THE Z-MOMENTUM EQUATION 
CALL PARTIALC1,W) 

CALL MOVLEVCDUDXC 1,1,1), GCI, 1,1), IJK) 

CALL PARTIAL ( 1 , G) 

DO 130 K=1 , LMAXP1 
DO 130 J=1 , JMAX 
DO 130 1=1, IMAX 

H3CI,J,K)=H3CI,J,K)+ExDUDXCI,J,K) 

130 CONTINUE 

CALL PARTIALC2, W) 

CALL MOVLEVCDUDXC 1,1,1),G(1,1,1),IJK) 

CALL PARTIAL C2, G) 
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DO 135 K=1 , LMAXP1 
DO 135 J=1 , JMAX 
DO 135 1=1 , IMAX 

H3CI, J,K)=H3CI,J,K)+EXDUDXCI,J,K) 

135 CONTINUE 
RETURN 
END 

XDECK SGS 

SUBROUTINE SGS 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

CXTHIS SUBROUTINE COMPUTES THE EDDY VISCOSITY AND THE SUBGRID SCALE X 

CXTERMS WHICH ARE ADDED TO THE RIGHT HAND SIDE OF THE GOVERNING MOMEN X 
CX-TUM EQUATIONS. THE EDDY VISCOSITY IS SET EQUAL TO ZERO AT THE WALL. X 
cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
COMMON/ADV/NTIME 

C0MM0N/SGTT/SGSTC65) , ETEDC65) ,U2STTC65) » V2STTC65 ) , W2STT < 65) 

1 » TSHGS, TSCNT 
COMMON/TINC/DT 
REAL MIU 

CQMMON/COUNT/IICONT 

COMMON/CON5T/C1 00 , C101,IJK,IJ> NHP1 > HALF 
COMMON/DATA9/IMAX, JMAX, LMAX, NHALFX, NHALFY 
COMMON/SCM2/LMAXP1,D1,D2,D3,D4,D5,D6 
COMMON/INNERC/CVINRC65) 

DIMENSION EDV0C65) , EDVI C65) 

XCALL A2 
XCALL A 9 
XCALL B2 
XCALL B3 
XCALL BA 
XCALL B5 
XCALL A4 
XCALL A7 
XCALL A6 
XCALL A5 

LMAXM1=LMAX-1 
IFCNTIME.NE.l) GO TO 5 
TSCNT=Q . 

TSHGS=0 . 

DO 2 K=1 , LMAXP1 
SGST(K)=0. 

ETEDCK)=0. 

U2STTCK)=0. 

V2STT ( K ) =0 . 

W2STT ( K ) =0 . 

2 CONTINUE 
5 CONTINUE 

LHPl=LMAX/2+l 

Cxxxxxx FIRST COMPUTE THE EDDY VISCOSITY, G. 

CALL PARTIAL( 1 , U) 

DO 10 K=3, LMAXM1 
DO 10 J =1, JMAX 
DO 10 1=1, IMAX 
G(I , J ,K)=DUDX(I , J,K)XX2 
10 CONTINUE 

CALL PARTIAL(2, V) 

DO 15 K=3, LMAXM1 
DO 15 J=l, JMAX 
DO 15 1=1, IMAX 

GCI , J,K)=G(I, J,K)+DUDX(I, J,K)XX2 
15 CONTINUE 

CALL PARTIALC3,W) 

DO 20 K=3 , LMAXM1 
DO 20 J=1 , JMAX 
DO 20 1=1, IMAX 

GCI, J,K)=GU,J,K)+DUDXCI, J,K)XX2 
20 CONTINUE 

CALL PARTIALCZ,U) 

CALL MOVLEV<DUDX( 1,1,1), PCI, 1,1), IJK) 

CALL PARTIALCl.V) 
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DO 25 K=3 , LMAXM1 
DO 25 J=1 > JMAX 
DO 25 I=1»IMAX 

GCI, J,K)=2.*GCI, J,K)+CDUDXCI,J,K)+PCI,J,K))X*2 
CONTINUE 

CALL PARTIALC2,W) 

CALL M0VLEVCDUDXC1,1> 1)»PC1,1>1)»IJK) 

CALL PARTIALC3,V) 

DO 30 K=3 > LMAXM1 
DO 30 J=1,JMAX 
DO 30 1=1, IMAX 

GCI,J,K)=GCI,J,K)+CDUDXCI,J,K)+PCI,J,K))X*2 

CONTINUE 

CALL PARTIALC1,W) 

CALL MGVLEVCDUDXC1,1,1),P(1,1,1),IJK) 

CALL PARTIALC3,U) 

BMAX=0. 

DO 35 K=3 , LMAXM1 
DO 35 J=1,JMAX 
DO 35 1=1, IMAX 

CCC=G(I,J,K)+CDUDXCI,J,K)+PCI,J,K))XX2 
H2(l, J,K)=CVCK)*SQRTCCCC) 

H1CI,J,K)=CVINRCK)*CCC 

CONTINUE 

CC=1 . / ( IMAX* JMAX) 

COMPUTE THE PLANAR AVERAGE OF INNER AND OUTER LAYER MODELS. 
DO 900 K=3, LMAXM1 
EDVGCK)=0 . 

EDVICK)=0. 

DO 910 J=1 , JMAX 
DO 910 1=1, IMAX 
EDV0(K)=EDV0CK)+H2CI, J»K) 

EDVICK)=EDVICK)+H1CI, J,K) 

CONTINUE 

EDVOC K) =EDV0CK)*CC 

EDVICK)=EDVI(K)*CC 

CONTINUE 

CR = 1 . 0 

MMM=0 

DO 915 K=3 , LHP1 

IFCEDVICK) . GT.EDVOCK)) MMM=2 

IFCMMM.EQ.2) GO TO 915 

IFCEDVICK). LT.EDVOCK)) KCROSl=K 

CONTINUE 

MMM = 0 

DO 920 K=LHP1 , LMAXM1 

KK=LMAXM1-K+LHP1 

IFC EDVI CKK) . GT . EDVOCKK) 1 MMM=2 

IFCMMM.EQ.2) GO TO 920 

IFC EDVICKK) -LT. EDVOCKK)) KCR0S2=KK 

CONTINUE 

PRINT 925 , KCROS1 , KCR0S2 

F0RMATC5X,* CROSS OVER POINTS OF INNER AND OUTER LAYER*, 215) 
PRINT 930 

FORMATC/,20X,* PLANE AVERAGE OF INNER LAYER MODEL *) 

PRINT 20 0 , C EDVI CK),K=3> LMAXM1 ) 

PRINT 935 

FORMAT C/ , 20X, X PLANE AVERAGE OF OUTER LAYER MODEL X) 

PRINT 200, C£DVOCK),K=3,LMAXMX) 

DO 940 K=3 , KCROS1 

DO 940 J=1 , JMAX 

DO 940 1=1, IMAX 

GCI,J,K)=H1CI,J,K)*CR 

CONTINUE 

KCR0S3=KCR0S1+1 

KCR0S4=KCR0S2-1 

DO 945 K=KCR0S3 , KCROS4 

DO 945 J=1 , UMAX 

DO 945 1=1, IMAX 

GCI,J,K)=H2CI»J»K) 
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945 CONTINUE 

DO 950 K=KCR0S2,LMAXM1 
DO 950 J=1 , JMAX 
DO 950 1=1 , IMAX 
GCI,J,K)=H1CI,J,K)XCR 
950 CONTINUE 

DO 4.0 J=l-> JMAX 
DO 40 1=1 , IMAX 
GCI,J,1)=0. 

GCI,J,2)=0. 

GCI,J,LMAX)=0. 

GCI,J,LMAXP1)=0. 

40 CONTINUE 

200 FORMATC IX, 1P9E14 . 5) 

CXXXXXX COMPUTE THE AVERAGE OF EDDY VISCOSITY IN X-Y PLANES 
DO 600 K=1 , LMAXP1 
MIUCK) =0 . 

DO 610 J=1,JMAX 
DO 610 1=1, IMAX 
MIUCiO=MIUCKH-GCI, J,K) 

610 CONTINUE 

MIUCK) =MIUCK)/ CIMAXX JMAX) 

600 CONTINUE 
PRINT 190 

190 FORMAT C10X,X AVERAGE EDDY VISCOSITY X) 

PRINT 200»C MIUCK), K~2> LMAX) 

CXXXXX COMPUTE THE VISCOUS INSTABILITY CRITERION. 

BMAX=0 . 

DO 400 K=3, LHP1 
KM1=K-1 

DO 400 J=1 , JMAX 
DO 400 1=1, IMAX 

VIS=CCZCK)-ZCKM1))XX2)/CABSCGCI, J,K)“MIUCK))) 

VIS=DT/VIS 

IFCVIS.LT. BMAX) GO TO 400 
BMAX=VIS 
IDUM2=I 
JDUM2=J 
KDUM2=K 
400 CONTINUE 
DMAX=0. 

DO 500 K=LHP1, LMAXM1 
KP1=K+1 

DO 500 J=1 , JMAX 
DO 500 1=1, IMAX 

VIS=( CZCKPD-ZCK) ) XX2 )/ CABS CGCI »J ,K) -MIUCK) ) ) 

VIS=DT/VIS 

IFCVIS.LT. DMAX) GO TO 500 
DMAX=VIS 
IDUM1=I 
JDUM1=J 
KDUM1=K 
500 CONTINUE 

PRINT 510, BMAX , I DUM1 , JDUM1 , KDUM1 , DMAX, IDUM2 , JDUM2 , KDUM2 
510 FORMAT C1X, X VIS INSTABILITY X,1P1E14.5,3I5,5X,1P1E14.5,3I5) 
CXXXXXEDDY VISCOSITY IS COMPUTED, NOW COMPUTE THE SUBGRID SCALE TERMS 
CALL PARTI ALC1 ,U) 

DO 60 K=l, LMAXP1 
DO 60 J=1 , JMAX 
DO 60 1=1, IMAX 

PCI,J,K)=2.XGCI,J,K)XDUDXCI,J,K) 

60 CONTINUE 

CALL PARTIALC1 ,P) 

DO 62 K=l, LMAXP1 
DO 62 J=1 , JMAX 
DO 62 1=1, IMAX 
HI C I , J » K) =DUDXCI , J » K) 

62 CONTINUE 

CALL PARTIALC2,U) 

CALL MOVLEVCDUDXC1 ,1,1)>PC1»1,1)>IJK) 


116 



original page is 
QE EGGS- QUALITX 


CALL PARTIALC1,V) 

DO 6 4 K=l, LMAXPI 
DO 64 J=l, JMAX 
DO 64 1=1 , IMAX 

PCI,J,K)=GCI, J,K)*CPCI, J,K)+DUDXCI,J,K)) 

64 CONTINUE 

CALL PARTIALC2,P) 

DO 66 K=1 » LMAXPI 
DO 66 J=1,JMAX 
DO 66 1=1, IMAX 

HI Cl , J,K)=H1CI, J,K)+DUDXCI, J,K) 

66 CONTINUE 

CALL PARTIALC3,U) 

CALL MO VL EVCDUDX (1,1,1), PCI, 1,1)»IJK) 

CALL PARTIALC1,W) 

DO 68 K=l, LMAXPI 
DO 68 J = 1 » JMAX 
DO 68 1=1, IMAX 

PCI,J,K)=PCI,J,K) +DUDX C I , J , K ) 

68 CONTINUE 

Cxxxxxx CALCULATE SOS CONTRIBUTIONS TO REYNOLDS STRESS AND INTENSITIES. 
CXXXXXX ALSO AVERAGE THEM IN TIME. 

TSHGS=TSHGS+1 
DO 92 K=l, LMAXPI 
SSUM ( K) =0 . 

DO 94 J=1 , JMAX 
DO 94 1=1, IMAX 

SSUMCK)=SSUMCK)+PCI,J,K)*GCI,J,K) 

94 CONTINUE 

SSUMCK) =-SSUMCK)/C IMAX* JMAX) 

SGST CK)=SGST(K)+$SUMCK) 

92 CONTINUE 

IFCNTIME.EQ.l) GO TO 360 
IFCIICONT.NE. 0) GO TO 350 
360 CONTINUE 

DO 98 K=l, LMAXPI 
EDYVICK)=0. 

DO 102 J=1 , JMAX 
DO 102 1=1, IMAX 
EDYVICK)=EDYVIOO+G(I,J,K)XX2 
102 CONTINUE 

EDYVI C K ) = EDYV I C K ) * FACTOR C K)/ C IMAX* JMAX) 

98 CONTINUE 

CALL PARTIAL C 1 , U ) 

DO 104 K=l, LMAXPI 
U2SCK)=0. 

DO 106 J=1 , JMAX 
DO 106 1=1, IMAX 

U2S CK)=U2S CK) + GC I, J ,K)XDUDXCI» J > K) 

106 CONTINUE 

U2SCK)=U2SCK)*2./CIMAX*JMAX) 

U2SCK)=EDYVICK)-U2SCK) 

104 CONTINUE ' 

CALL PARTIALC2,V) 

DO 108 K=l, LMAXPI 
V2SCK) =0 . 

DO 110 J=1 , JMAX 
DO 110 1=1, IMAX 

V2SCK)=V2SCK)+GCI,J,K)*DUDXCI,J,K) 

110 CONTINUE 

V2S CK)=V2S CK)X2 ./C IMAXXJMAX) 

V2SCK)=EDYVICK)-V2SCK) 

108 CONTINUE 

CALL PARTIALC3,W) 

DO 112 K=l, LMAXPI 
W2SCK)=0 . 

DO 114 J =1 , JMAX 
DO 114 1=1, IMAX 

W2SCK)=W2S(K)+GCI, J,K)*DUDXCI,J,K) 

114 CONTINUE 


117 



W2S C K > =W2S C K) X2 . / C IMAX* JMAX) 
W2SCK)=EDYVICK)-W2S<K) 

112 CONTINUE 

TSCNT=TSCNT+1 

DO 220 K=3,LMAXM1 

ETED(K)=ETEDCK)+EDYVI^K) 

U2ST-T-(-K)=U2STTCK>+U2SaC) 

V2STT(K)=V2STT (K)+V2S(K) 

W2STTCK ) =W25TT ( K) +W2S(K) 

220 CONTINUE 
350 CONTINUE 

CALL PARTIALC3,P) 

DO 70 K=1 , LMAXP1 
DO 70 J=1,JMAX 
DO 70 1=1, IMAX 

Hia,J,K)=Hia,J,K)+(G(I,J,K>-MIU(K))*DUDX(I,J,K) 

70 CONTINUE 

CALL PARTIAL(3> G) 

DO 71 K=1,LMAXP1 
DO 71 J =1 , UMAX 
DO 71 1=1, IMAX 

HlCI,J,K)=Hl<I,J,K)+DUDXCI,J,K)XPa,J,K) 

71 CONTINUE 

CALL PARTIAL ( 1 , W) 

CALL MOVLEV(DUDX(l»l,l).»Ptl»l»l),IJK) 

CALL PARTIALC3,P) 

DO 715 K=1 > LMAXP1 
DO 715 J=1,JMAX 
DO 715 1=1, IMAX 

Hl(I,J,K)=HlCI,J,lO+MIUCiOXDUDXa,J,iC> 

715 CONTINUE 

CXXXXX Y-MOMENTUM EQUATION. 

CALL PARTIALC1, V) 

CALL M0VLEV(DUDXC1,1,1),P(1,1,1),IJK) 

CALL PARTIAL(2,U) 

DO 72 K = 1 » LMAXP1 
DO 72 J=1,JMAX 
DO 72 1=1, IMAX 

P(I,J,K)=G( I, J,K)X(P(I, J,K)+DUDX(I> J,K)) 

72 CONTINUE 

CALL PARTIALC1, P) 

CALL M0VLEV(DUDXCl,l,l),H2a,l,l),IJK) 

CALL PARTIAL(2» V) 

DO 74 K = l» LMAXP1 
DO 74 J =1 , JMAX 
DO 74 1=1, IMAX 

P(I,J,K)=2.XGCI,J,K)XDUDXa,J,K) 

74 CONTINUE 

CALL PARTIALC2,P) 

DO 76 K=l, LMAXP1 
DO 76 J =1 , JMAX 
DO 76 1=1, IMAX 

H2(I,J,K)=H2(I,J,K)+DUDXCI,J,K) 

76 CONTINUE 

CALL PARTIALC3, V) 

CALL MOVLEVCDUDXC1,1,1),PC1,1,1),IJIO 
CALL PARTIAL(2,W) 

DO 78 K=l, LMAXP1 
DO 78 J=l» JMAX 
DO 78 1=1, IMAX 

P(I,J,K>=P<I,J,K)+DUDXa,J,K> 

78 CONTINUE 

CALL PARTIAL(3,P) 

DO 80 K=1 , LMAXP1 
DO 80 J=l, JMAX 
DO 80 1=1, IMAX 

H2(I»J,K)=H2(I,J , K)+(G(I» J , K)-MIU(K) ) XDUDXC I , J , K) 
80 CONTINUE 

CALL PARTIAL (3 , G) 

DO 81 K=l, LMAXP 1 
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DO 81 J=1,JMAX 
DO 81 I=1,IMAX 

H2a,J,K3=H2a,J,K3+DUDXCI,J,K)XP(I,J,K) 

81 CONTINUE 

CALL PARTI AL(2,W3 

CALL M0.VLEVCDUDXC1,1,13,P(1,1,13,IJK3 
CALL PARTIAL(3,P3 
DO 815 K=1 » LMAXPI 
DO 815 J=1,JMAX 
DO 815 I = 1,-IMAX 

H2a,J,K)=H2(I,J,K)+MIU(K)XDUDXa, J,K) 

815 CONTINUE 

CXXXXXX Z-MOMENTUM EQUATION. 

CALL PARTIALC1,U3 

CALL MOVLEVCDUDXCl, 1,1 ),P(1,1»1) »I JK) 

CALL PARTIAL(3,U) 

DO 82 K=1 , LMAXPI 
DO 82 J=1,JMAX 
DO 82 1=1, IMAX 

P(I, J,K3=GCI, J,K3XCPCI,J,K3+DUDX(I,J,K>) 

82 CONTINUE 

CALL PARTIALC1,P) 

CALL MOVLEVCDUDXC 1 ,l,l),H3(l,l,l),IJf() 

CALL PARTI AL(2,W3 

CALL M0VLEVCDUDXC1,1,13>PC1,1,1),IJK) 

CALL PARTIAL(3,V3 
DO 84 K=1 , LMAXPI 
DO 84 J=l, JMAX 
DO 84 I=1,IMAX 

Pa,J,K3=G(I,J,K3*CPCI,J,K3 + DUDXCI,J»K)3 
84 CONTINUE 

CALL PARTIAL(2,P3 
DO 86 K=l, LMAXPI 
DO 86 J=l, JMAX 
DO 86 1 = 1, MAX 

H3U, J,K3 = H3(I, J,K3+DUDX(I,J,K3 
86 CONTINUE 

CALL PARTIAL C3,W3 
DO 88 K=l, LMAXPI 
DO 88 J=l, JMAX 
DO 88 1=1 , IMAX 
PCI, J,K)=2.*DUDX(I,J,IO 
88 CONTINUE 

CALL PARTIALC3,P) 

DO 90 K=l, LMAXPI 
DO 90 J=1 , JMAX 
DO 90 1=1, IMAX 

H3(I,J,K)=H3a,J,K) + CGCI,J,K)“MIUCK))XDUDXCI,J,K) 

90 CONTINUE 

CALL PARTIALC3,G3 
DO 91 K=l, LMAXPI 
DO 91 J=1 , JMAX 
DO 91 1=1, IMAX 

H3CI, J,K)=H3CI, J,K3+DUDXCI,J,K3XPCI, J,K) 

91 CONTINUE 

DO 100 J = 1 , JMAX 
DO 100 1=1, IMAX 
Hia,J,l) = 0. 

H1CI,J,2)=0. 

H2CI , J, 13=0 . 

H2CI,J,23=0. 

H3(I,J,13=0. 

H3(I,J,23=0. 

HI ( I , J , LMAX3=0 . 

H2 (I , J , LMAX3 =0 . 

H3 ( I , J , LMAX3 =0 . 

HI C I , J , LMAXPI 3 = 0 . 

H2(I,J, LMAXPI 3=0. 

H3 ( I , J , LMAXPI 3 = 0 . 

100 CONTINUE 
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RETURN 

END 

XDECK FILTER 

SUBROUTINE FILTER(HR) 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXX 

CX THIS SUBROUTINE FILTERS A THREE DIMENSIONAL ARRAY IN X AND Y DIRECTIONS 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

COMMON/DATA9/IMAX, JMAX, LMAX,NHALFX, NHALFY 
COMMON/ SCM2/LMAXP1 , D1 , D2, D3 , DA , D5 , D6 
COMMON/CONST/CIOO, Cl 01 , I JK, I J »NHP1 , HALF 
XCALL CI1 
XCALL A1 
XCALL B1 

CC=1./CIMAXKJMAX) 

DO 20 L=1,LMAXP1 

CALL MOV LEV (HR(1»1»L)»FRC1,1)>IJ) 

CALL FFTXC1.05 
CALL FFTYC1. 0,1.0) 

DO 30 1=1 > IMAX 
DO 30 J=1,JMAX 

FR(I, J)=FR(I, J)XFILTX(I)XFILTYCJ) 

Fia,J)=FItI,J)XFILTXCI)*FILTY(J) 

30 CONTINUE 

CALL FFTY (-1 . 0 >CC) 

CALL FFTXC-I.O) 

CALL M0VLEVCFRC1,1),HRC1,1,L),IJ) 

20 CONTINUE 
RETURN 
END 

XDECK STAT 

SUBROUTINE STAT 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

CX THIS SUBROUTINE COMPUTES THE STATISTICS OF THE FLOW FOR OUTPUT. X 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

COMMON/CONST/C1QO,C101,IJK,IJ,NHP1,HALF 
COMMON/SCM2/LMAXP1 , D1 , D2 , D3 , DA , D5 , D6 
COMMON/SCMA/CI,CJ,CK,CJK,CIK,CIJ 
COMMON/DATA9/IMAX, JMAX, LMAX, NHALFX, NHALFY 
XCALL BA 
XCALL B6 
XCALL A6 
XCALL A 9 

PRINT 2000 
2000 FORMATC1H1) 

PRINT 1100 

1100 FORMAT ClX t x UAVG IN X-Yx, AX,XVAVG IN X~YX, 3X,XWAVG IN X-YX,1X,X 
1U2AVG IN XYX,3X,XV2AVG IN XYX,3X,*W2AVG IN XY*,3X,XQ2AVG IN XYX 
1 » 3X» XTURB SHEARX,7X,XZX) 

UTOT=0. 

VTOT=0. 

WTOT=0. 

U2TOT=0. 

V2TOT=0. 

W2T0T=0. 

QTOT=0. 

PAI=ACOS(-l . ) 

DO 100 K=1 » LMAXP1 
USUM=0. 

VSUM=0. 

WSUM=0. 

DO 110 J=l> JMAX 
DO 110 1=1, IMAX 
USUM=USUM+U ( I , J , K) 

VSUM=VSUM+V(I, J,K) 

WSUM=WSUM+W( I » J > K) 

110 CONTINUE 

USUM=USUMXCIJ 
VSUM=VSUMXCIJ 
WSUM=WSUMXCIJ 
SHEAR=0 . 
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U2SUM=0. 

V2SUM=0. 

W2SUM=0. 

DO 120 J =1 > JMAX 
DO 120 I=1,IMAX 

U2SUM=U2SUM+(UCI, J,K)-USUM)XX2 
V2SUM=V2SUM+(VCI,J,K)-VSUM)XX2 
W2SUM=W2SUM+(W(I> J »K)-WSUM)XX2 
SHEAR=SHEAR+CUCI, J,K)-USUM)*(W(I, J,K)-WSUM) 

120 CONTINUE 

Q=(U2SUM+V2SUM+U2SUM)XCIJX.5 

U2SUM=SQRTCU2SUMXCIJ) 

V2SUM=SQRTCV2SUMXCIJ) 

W2SUM=SQRT(W2SUMXCIJ> 

SHEAR=SHEARXCIJ 

PRINT 1000,USUM,VSUM,WSUM,U2SUM,V2SUM,W2SUM,Q,SHEAR,ZCK) 
U2ST(IO=SQRT<U2SUM*X2+U2S<K)) 

V2ST(K)=SQRT(V2SUMXX2+V2S(Kn 

W2STCK)=SQRTCW2SUMXX2+W2S(K)) 

UWT(K) =SHEAR+5SUM(K) 

UT0T=UT0T+USUM 
VT0T=VTOT+VSUM 
WT0T=WT0T+WSUM 
U2T0T=U2T0T+U2SUM 
V2TOT=V2TOT+V2SUM 
W2T0T=W2T0T+W2SUM 
QTOT=QTOT+Q 
100 CONTINUE 

UTOT=UTOTXCK 

VTOT=VTOT*CK 

WTOT=WTOTXCK 

U2T0T=U2T0TXCK 

V2T0T=V2T0TXCK 

W2T0T=W2T0TXCK 

QTOT=QTOTXCK 

PRINT 1200 

1200 FORMATC///,lX,X UTOT IN X-Y VTOT IN X-Y WTOT IN X-Y U2T0T 
1 IN X-Y V2T0T IN X-Y W2T0T IN X-Y TURB ENERGY *5 
PRINT 1000, UTOT , VTOT, WTOT, U2T0T , V2T0T , W2T0T , QTOT 

1000 FORMAT (1P9E14.5) 

PRINT 200 

200 FORMATC//, 5X, X INSTANTENEOUS UX) 

PRINT 210>(U(8»8»K),K = 1» LMAXP1 ) 

210 FORMAT (1X»1P9E14.5) 

PRINT 300 

300 FORMATC////, 30X, X SGS CONTRIBUTIONS ADDEDX) 

PRINT 310 

310 FORMAT C IX, X SGS ENERGYX, AX,x T0TALU2S x, 5 X,x TOTAL V2S X,3X,X TOT 
1ALW2S x,3X,XT0TAl SHEARX,3X,x PLANEX) 

LMAXM1=LMAX-1 
DO 320 K=3, LMAXM1 

PRINT 330 , EDYVI ( K) , U2STCK) , V2STGO , W2STCK) ,UWTCK) , K 
320 CONTINUE 

330 FORMAT ( IX >1P5 El 4.5,16) 

RETURN 

END 

XDECK TRANS 

SUBROUTINE TRANS 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

cx THIS SUBROUTINE COMPUTES THE VARIOUS TRANSFORMATION QUANTITIES X 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

C0MM0N/DATA9/IMAX, JMAX, LMAX,NHALFX,NHALFY 
C0MM0N/LENGTH/LSCALEC65) 

REAL LSCALE 

C0MM0N/INNERC/CVINRC65) 

C0MM0N/SCM3/DELTA1 , DELTA2, RE, E 

COMMON/RANGE/LMAXM1 , LMAXM2, LMAXM3 , LMAXM4 , LMAXM5 
COMMON/TINC/DT 

C0MMON/BC/CE1 , CE2, CE3»CE4, CE5,CE6 
XCALL B2 
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XCALL B3 
XCALL A 9 
XCALL C7 
XCALL B7 
XCALL B8 

COMMON/PENTA2/XI,QI,GI,YI,QJ,GJ,XN,QIN,GIN,YN,.QJN,GJN,Q2,Q3,- 
1RC1,RC2,_RP1,RR2,RP-3,RP4 
C0MM0N/ZER0/C3 , CA 
COMMON/IDENTN/CODE 
LMAXM1=LP1AX-1 
LMAXM2=LMAX-2 
LMAXM3=LMAX-3 
LMAXMA=LMAX-A 
LMAXM5=LMAX-5 
LMAXP1=LMAX+1 
LHPl=LMAX/2+l 

CXXXXX MESH STRECHING TRANSFORMATION 
P=0. 983A6 

TANIP=0.5XALOG((1. +?)/(!. - P)) 

PINV=1./P 

P2=PX*2 

DO 5 J = 1 > LMAXP1 

ZETAU)=-l. + 2.X(J-2)/CLMAX-2) * 

DUMl-ZETAC J )XTANIP 
ZCJ)=PINVXTANH(DUM1) 

RL(J)=C2.XP2/TANIP)XCCC0SH(DUM1))X*3)X(SINHCDUM1)) 

RMCJ)=P2X((C0SHCDUMl>)xx4>/aANIPxx2) 

5 CONTINUE 

DELTA3=ZETA(2)-ZETA(1 ) 

E2=RL(2)/(2.XDELTA3XRE> 

F2=RMC2)/<CDELTA3XX2)XRE) 

EN=RLCLMAX)/(2.XDELTA3XRE) 

FN=RMCLMAX)/C (DELTA3X*2)XRE) 

R2=( F2+E2)/ ( F2-E2 ) 

RN=( FN-EN)/( FN+EN) 

RR2=1 ./( E2-F2 ) 

RRN=-1 ./(EN+FN) 

PRINT 20 

20 F0RMATC6X,XZETAX,12X,*ZX,14X,*RL*,1AX,XRMX) 

DO 30 K=1 > LMAXP1 

PRINT AO,ZETACK),Z(K),RLCK),RM<K) 

30 CONTINUE 

AO FORMAT (1X,1P4E15.7) 

PRINT 5Q,E2,F2,EN,FN,R2,RN,DELTA3 
50 F0RMAT(1X,//,1P7E1A.5) 

CC=0 . 2 

CXXXXX COMPUTE THE LENGTH SCALE FOR THE SGS MODEL 
VONK=O.A 

DFILT1=2 .XDELTA1 
DFILT2=2.XDELTA2 
P0WER=l./3. 

DO 300 K=3,LHP1 
KM1=K-1 

DW=(ZCK)-Z(2))XV0NK 

GRID=ZCK)~Z(KM1) 

LSCALE(K) =( AMIN1 (DW> 0 . 1, DFILT1) )X(AMIN1(DW, 0 . 1»DFILT2) )X(AMIN1 
1(DW, 0 . 1 ,GRID) ) 

LSCALE(K) ~LSCALECK)XXPOWER 
300 CONTINUE 

DO 310 K=LHP1 » LMAXM1 
KK=LMAXM1-K+LHP1 
KKP1 =KK+1 

DW=CZCLMAX)-Z(KK))XVONK 

GRID=ZCKKP1)-Z(KK) 

LSCALE(KK) = (AMIN1(DM,0.1,DFILT1))X(AMINKDM,0.1,DFILT2))X<AMIN1 
1CDW, 0.1, GRID)) 

LSCALE(KK)=LSCALECKK)XXP0WER 
310 CONTINUE 

CINER=CCCXX2)/CV0NKX27.) 

DO 320 K=3 , LMAXM1 
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CVCK) = CCCXLSCALECK) )XX2 
CVINRCK)=CINERXREXCLSCALECK))XX4 
320 CONTINUE 
PRINT 330 

330 FORMAT C// ,20X,X COEFFICIENT OF INNER SGSX) 

PRINT 120, CCVINRCK),K=3, LMAXM1) 

PRINT 340 i 

340 FQRMATC//,20X,x SUBGRID LENGTH SCALEX) 

PRINT 120,CLSCALECK),K=3» LMAXM1 ) 

PRINT 110 

110 FORMATC 20 X, x COEFFICIENT OF SGS x) 

PRINT 120, CCVCK)»K = 3> LMAXM1) 

120 FORMAT C1X,1P9E14.5) 

FAC=226X(CCXX2)/3. 

FACTORC1 ) =0 . 

FACTORC2) =0 . 

FACTORC LMAX) = 0 . 

FACTORC LMAXP1 ) =0 . 

DO 100 K = 3 , LMAXM1 
FACTORCK) =FAC/CVCK) 

100 CONTINUE 

DO 12 J=2 , LMAX 
H1=ZCJ)-ZCJ-1) 

H2=ZCJ+1)-ZCJ> 

CXXXXX ARRAYS FOR FINITE DIFFERENCE IN Z-DIRECTION 
AP C J ) =-l . /CZC J+l )-ZC J-l) ) 

BPCJ)=0. 

CPCJ)=-APCJ) 

Cxxxxxxxxxxx DEFINE THE COEFFICIENTS FOR SECOND DERIVATIVE IN Z DIREC 
AP2CJ)=2./(H1*(H1+H2)> 

BP2CJ)=-2./(HlXH2) 

CP2CJ)=2./CH2XCH1+H2)) 

PRINT 80,APCJ),BPCJ),CPCJ)> AP2U) , BP2C J ) ,CP2C J ) 

12 CONTINUE 

Cxxx CONSTANTS FOR THE BLOCK TRI-DIAGONAL MATRIX IN THE MAIN PROGRAM 
T=0.5X(Z(3)-ZC2)) 

CEl=l.-AP(3)xTXEXDTX0.5X(CP2(2)-AP2(2)XCP(2)/AP(2n/a.+TXAPC3)) 

CE2=BPC3)+AP(3)X(I.-TXBPC3))/(1.+TXAPC3)) 

CE3=CP(3)-AP(3)XTXCP(3)/(1.+TXAPC3)) 

T = 0.5X(Z( LMAX)-ZC LMAXM1)) 

CE4=AP( LMAXM1 ) + CP( LMAXM1 JXTXAPC LMAXM1 )/( 1 TXCPC LMAXM1 ) ) 

CE5=BP C LMAXM1) +CPC LMAXM1 )X( 1 . +TXBPC LMAXM1) )/C 1 . -TXCPC LMAXM1 ) ) 

CE6=1.+CPCLMAXM1)XTXEXDTXG.5X(AP2CLMAX)-CP2UMAX>*APCLMAX) 

1/CPC LMAX) )/(l .-TXCPC LMAXM1) ) 

T=0.5XCZC3)-ZC2) ) 

C3=C1.-TXBPC3))/CPC3) 

C4=CTXCPC3)/(1.-TXBPC3))) 

Q=1./C1.+TXAPC3)) 

XI=-TXQ 

QI = C 1 . -TXBP C 3) ) XQ 

GI=-TXCPC3)XQ 

YI=C1.+BPC2)XTXQ)/APC2) 

QJ=CBPC2)XCTXBPC3)-1. )XQ-CPC2))/APC2) 

GJ=BPC2)XTXCPC3)XQ/APC2) 

T=0 . 5X (Z( LMAX)-ZCLMAXM1 ) ) 

Q=1 ./(l . -TXCPC LMAXM1 ) ) 

XN=TXQ 

QIN=C1.+TXBPCLMAXM1))XQ 

GIN=TXAPCLMAXM1)XQ 

YN=C 1. -BP C LMAX) XTXQ)/CPC LMAX) 

QJN=-CBPCLMAX)xCl.+TxBPCLMAXMl))xQ+AP(LMAX))/CPCLMAX) 

GJN=-TXAP C LMAXM1 ) XBP C LMAX)XQ/CP C LMAX) 

80 FORMAT C IX, 1P3E15.7 ,5X,1P3E15.7) 

90 FORMAT C IX, IPS El 5. 7) 

RETURN 

END 

XDECK VISCOS 

SUBROUTINE VISCOSCU) 

Cxxxxx THIS SUBROUTINE COMPUTES THE SECOND DERIVATIVE OF U IN THE Z-DIRECTION 
COMMON/DA T A 9/IMAX, JMAX, LMAX, NHALFX, NHALFY 
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XCALL A2 
XCALL B7 
XCALL B 9 
XCALL A 9 

LMAXP1=LMAX+1 
DELTA3=2./(LMAX-2. ) 

DO 20 J = 1-,JMAX 
DO 20 1=1, IMAX 
DUDXCI,J,1)=0. 

DUDXC I , J , LMAXP1 ) =0 . 

20 CONTINUE 

DO 30 K=2,LMAX 
DO 30 J=1 , JMAX 
DO 30 1=1 , IMAX 
KP1-K+1 
KM1=K-1 

DUDXCI,J,K)=AP2(iOxua,J,KMl)+BP2CK>xua,J,K)+CP2CK)XU(I,J,KPl) 

30 CONTINUE 
RETURN 
END 

#DECK EXTERN 

SUBROUTINE EXTERNC LI , L2, R,RR) 

C XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

C X THIS SUBROUTINE FIXES THE EXTERNAL VALUES OF THE U AND V AND W X 

C X NOTE THAT THE EXTERNAL VALUES OF U AND V WILL NOT ENTER INTO THE X 

C X COMPUTATION. AND THEY ARE UNNECESSARY X 

C xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
COMMON/CONST/CIO 0 , CIO! , I JK , I J , NHP1 , HALF 
COMMON/S CM3/DELT A 1 >DELTA2» RE, E 
COMMON/DATA 9/IMAX, JMAX»LMAX,NHALFX,NHALFY 
XCALL A1 
XCALL C3 
XCALL A6 
XCALL C7 

LMAXP1=LMAX+I 
LMAXM1=LMAX-1 
DO 90 J = 1 » JMAX 
DO 90 1=1, IMAX 

WCI,J,1)=-CPC2)XUCI,J,3)/AP(2) 

90 W ( I , J , LMAXP 1 ) =- AP ( LMAX) XWC I , J , LMAXM1 )/CP C LMAX) 

DO 97 J=1 > JMAX 
DO 97 1=1, IMAX 
U(I,J,1)=0. 

VCI,J,1)=0. 

U(I, J,LMAXP1)=0. 

VC I , J , LMAXP1) =0 . 

97 CONTINUE 
95 CONTINUE 
RETURN 
END 

XDECK MTDAG 

SUBROUTINE MTDAGC AM, A, AP , V, N, K) 

C SOLVES COUPLED TRI-DIAGONAL ALGEBRAIC EQUATIONS A. 

C AM(I,J,L)X(J,L-1)+ACI»J,L)X(J, L )+AP(I, J » L )X( J , L+1)=YCI,L) 5. 

C (SUM OVER J IN EACH EQUATION) 6. 

C ( I , J , L ) I IS EQUATION TYPE, J IS VARIABLE TYPE, L IS NODE 7. 

C AT CALL VCI,L)=X(I,L) Y(J,L) IS RETURNED IN V(J,L) 8. 

C THE AM,A,AP ARRAYS ARE RETURNED AS GARBAGE 9. 

REAL AM(N,N,K) ,A(N,N,K),AP(N,N,K),V(N,K) 10. 

COMMON/SING/IMR, JMR, IMI , JMI 

C ELIMINATE TO OBTAIN A SEQUENTIALLY SOLVABLE FORM 11. 

DO 20 LX=1,K 12. 

L=K-LX+-1 13. 

LM=L-1 14. 

DO 18 J=1,N 15. 

C=ACJ, J,L) 16. 

IF CC.EQ.O.) GO TO 80 17. 

C ELIMINATE X(J,L) FROM ALL EQUATIONS OTHER THAN ITS OWN 18. 

DO 16 1=1, N 19. 

C ELIMINATE XCJ,L) FROM THE EQUATION FOR THE NODE L-l 20. 
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IF (L.EQ.l) GO TO 12 21. 

F=AP(I,J,LM> 22. 

IF (F.EQ.0.0) GO TO 12 23. 

F=F/C 24. 

DO 6 J1=1,N 25. 

A(I,J1,LM)=A(I,J1,LM)-FXAM(J,J1,L) 26. 

6 APCI,J1,LM)=AP(I,J1,LM)-F*A(J,J1,L> 27. 

V(I,LM)=V(I,LM)-F*VCJ,L) 28. 

C ELIMINATE XCJ,L) FROM OTHER EQUATIONS AT THIS NODE L 29. 

12 IF (I.EQ.J) GO TO 16 30. 

F=ACI,J,L) 31. 

IF (F.EQ.0.0) GO TO 16 32, 

F=F/C • 33. 

DO 14 J1=1,N 34. 

ACI,J1,L)=ACI,J1,L)-F*ACJ,J1,L> 35, 

IF (L.EQ.l) GO TO 14 36. 

AM(I,J1,L)=AM(I, J1,L)-F*AM(J,J1,L) 37. 

14 CONTINUE 38. 

V(I,L)=V(I,L)-F*VCJ,L> 39. 

16 CONTINUE 40. 

18 CONTINUE 41. 

20 CONTINUE 42. 

C CARRY OUT THE BACK SOLUTION 43. 

DO 30 L=1,K 44. 

LM=L-1 45. 

DO 28 1=1, N 46. 

C=A(I,I,L) 47. 

IF (C.EQ.0.0) GO TO 80 48. 

F=V(I,L) 49. 

IF (L.EQ.l) GO TO 28 50. 

DO 24 J1=1,N 51. 

24 F=F-AM(I, J1,L)*V(JI,LM) 52. 

28 V(I,L)=F/C 53. 

30 CONTINUE 54. 

RETURN 55. 


80 PRINT 90 

PRINT 10,IMR,JMR,IMI, JMI 
10 F0RMAT(4X,4I5> 

RETURN 

90 FORMAT (/// » 10X» * MTDAG MATRIX IS SINGULAR X) 

END 59. 

XDECK DIVG 

SUBROUTINE DIVG 

C THIS SUBROUTINE COMPUTES THE DIVERGENCE OF VELOCITY FIELD 
COMMON/DATA 9/IMAX, JMAX, LMAX,NHALFX,NHALFY 
COMMON/CONST/CIO 0 , C101,IJK,IJ ,NHP1 , HALF 
XCALL A2 
XCALL A6 
XCALL A5 

CALL PARTIAL(1,U) 

CALL MOVL EVCDUDXC 1 ,1,1),G(1,1,1)»IJK) 

CALL PARTIAL(2, V) 

DO 10 K=2,LMAX 
DO 10 J=1 » JMAX 
DO 10 1=1 , IMAX 

G(I,J,K)=G(I,J , K)+DUDX(I, J , K) 

10 CONTINUE 

CALL PARTIAL(3,W) 

DO 20 K=2, LMAX 
DO 20 J=1 , JMAX 
DO 20 1=1, IMAX 

G(I, J,K)=G(I,J,K)+DUDXCI,J,K) 

20 CONTINUE 
BMAX=0. 

DO 30 K=2, LMAX 
DO 30 J = 1 , JMAX 
DO 30 1=1, IMAX 

IF(ABS(G(I,J,K)> .GT.BMAX) BMAX=ABS(G(I, J , K) ) 

30 CONTINUE 

PRINT 40 , BMAX 
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40 F0RMAT(2X,X MAX DIVERGENCE^*, 1P1E15. 7) 

RETURN 

END 

XDECK COURANT 

SUBROUTINE COURANT (DT,NTIME, TEND) 

CXXXXX THIS SUBROUTINE MONITORS THE COURANT NUMBER 
XCALL A 9 
XCALL A 6 

C0MM0N/SCM3/DELTA1 ,DELTA2, RE, E 

C0MM0N/DATA9/IMAX, JMAX, LMAX, NHALFX, NHALFY 

LMAXM1=LMAX-1 

LHPl=LMAX/2+l 

BMAX=0. 

DO 51 K=3 , LHP1 
KM1=K-1 

DO 51 J=1 , JMAX 
DO 51 1=1 , XMAX 

CMAXl=ABSCWa,J,K))/(ZCK)-Z(KMl))+ABSCU(I,J,K)/DELTAl)+ABS(V(I,J, 
IK) )/DELTA2 

IFCCMAX1 . LT.BMAX) GO TO 51 
BMAX=CMAX1 
IDUM1=I 
JDUM1=J 
KDUM1=K 
51 CONTINUE 
DMAX=0. 

DO 56 K=LHP1,LMAXM1 
KP1=K+1 

DO 56 J = 1 > JMAX 
DO 56 I“1 , IMAX 

CMAX2-ABS(W(I» J»K) )/(Z(KPl)-Z(K))+ABS(U(I» J»K))/DELTA1+ABS( 
1V(I,J,K))/DELTA2 
IF(CMAX2. LT.DMAX) GO TO 56 
DMAX=CMAX2 
IDUM2=I 
JDUM2=J 
KDUM2=K 
56 CONTINUE 
BMAX=BMAXXDT 

dmax=dmaxxdt 

PRINT 61 , BMAX, IDUM1 , JDUM1 , KDUM1 ,DMAX, I DUM2 , JDUM2 » KDUM2 
61 F0RMATC2X,* COURRANT X,1P1E14.5,3I5,1P1E14.5,3I5> 

I FC BMAX . GT . 0 . 35 . OR . DMAX . GT . 0 . 35 ) NTIME=TEND 

RETURN 

END 

XDECK LTAVG 

SUBROUTINE LTAVG 

CXXXXX THIS SUBROUTINE COMPUTES THE RUNNING TIME AVERAGE OF VARIOUS 
CXXXXX STATISTICAL QUANTITIES. 

CQMM0N/SCM4/CI , CJ, CK, CJK, CIK, Cl J 
COMMON/DATA 9/IMAX, JMAX, LMAX, NHALFX, NHALFY 
XCALL A2 
XCALL A3 
XCALL A4 
XCALL A5 
XCALL A6 

COMMON/RANGE/LMAXM1 , LMAXM2, LMAXM3, LMAXM4, LMAXM5 
C0MM0N/SCM2/ LMAXP1 , D1 , D2 , D9 , D4 , D5, D6 

COMMON/LT A1/USUM(65) , UTSUMC65) , STSUM(65) »U2SMT{65) , V2SMT(65) 

1 >W2SMT(65) , PVT (65) , PUT(65) ,PUNST(65) ,PVNST(65) ,PMNST(65) , PWT(65) 
2,TC0NT 

C0MM0N/LTA2/PDUTC65) ,PDVT(65), PDWT(65) ,PDUNT(65) , PDVNT(65) ,PDWNT 
1(65) 

COMMON/ADV/NTIME 
IFCNTIME.NE.l) GO TO 5 
TC0NT=0. 

DO 2 K=3, LMAXM1 
UTSUM(K) =0 . 

U2SMT(K)=0. 

V2SMT(K)=0. 
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W2SMT C K ) =0 . 

STSUMCK ) =0 . 

PUT CK) =0 . 

PVTCK)=0 . 

PWT (K)-0 . 

PUNSTCK)=0. 

PVNST CK)=0 . 

PWNSTCK)=0 . 

PDUTCK) =0 . 

PDVT CK) =0 . 

PDWTCK) =0 . 

PDUNT CK)=0 . 

PDVNTCK)=0. 

PDWNTCK)=0. 

2 CONTINUE 
5 CONTINUE 
TCQNT=TCONT+l 
DO 10 K=3 > LMAXM1 
USUM(K)=0 . 

DO 15 J=l> UMAX 
DO 15 1=1, IMAX 
USUMCK)=USUMCK)+UCI,J,K> 

15 CONTINUE 

USUM(K)=USUMCK)XCIJ 

UTSUMCK)=UTSUMCK)+USUMCK) 

10 CONTINUE 

DO 20 K=3, LMAXM1 
U2SUM=0. 

V2SUM=0. 

W2SUM-0 . 

SSUM=0. 

DO 25 J=1,JMAX 
DO 25 1=1, IMAX 

U2SUM=U2SUM+ ( U (I , J , K) -USUMCK) )XX2 
V2SUM=V2SUM+VCI,J,K)XX2 
W2SUM=W2SUM+WCI, J,K)X*2 
SSUM=SSUM+WCI,J,K)XCUCI,J,K)-USUMCK>) 

25 CONTINUE 

U2SUM=U2SUM*CIJ 
V2SUM=V2SUM*CIJ 
W2SUM=W2$UMXCIJ 
SSUM=SSUMXCIJ 
U2SMTCK) =U2SMT (K1+U2SUM 
V2SMTCK)=V2SMTCK)+V2SUM 
W2SMTCK) =W2SMT CK) +W2SUM 
STSUMCK) =$TSUMCK)+SSUM 
20 CONTINUE 

DO 30 K“1 , LHAXP1 
DO 30 J=1,JMAX 
DO 30 1=1, IMAX 

PCI, J,K)=CUCI,J,K)x*2+VCI,J,K)xx2+WCI,J,K)XX2)/2. 
30 CONTINUE 

CALL FILTER(P) 

DO 35 K=1 , LMAXP1 
DO 35 J=1,JMAX 
DO 35 1=1, IMAX 

DIVCCI,J,K)=GCI,J,K)-PCI,J,K) 

35 CONTINUE 

CALL PARTI ALC1,DIVC) 

DO 40 K=3, LMAXM1 
PU=0. 


DO 4 5 J=l, JMAX 
DO 45 1=1, IMAX 
PU=PU+DUDXCI,J,K)XUCI, J,K) 
45 CONTINUE 
PU=PUXCIJ 
PUTCK)=PUTCK)+PU 
40 CONTINUE 

CALL PARTIAL Cl , G) 

DO 50 K=3, LMAXM1 
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PUNS=0. 

DO 55 J=l, JMAX 
DO 55 1=1, IMAX 

PUNS=PUNS+DUDXCI,J,K)XU(I,J,K) 
55 CONTINUE 

PUNS=PUNS*CIJ 
PUNSTCKO =PUNST(K)+PUNS 
50 CONTINUE 

CALL PARTIAL(2,DIVC) 

DO 60 K=3,LMAXM1 
PV=0. 

DO 65 J=1 > JMAX 
DO 65 1=1 , IMAX 
PV=PV+DUDX(I, J,K)*V(I> J>K) 

65 CONTINUE 
PV=PV*CIJ 
PVT(K) =PVT(K)+PV 
60 CONTINUE 

CALL PARTIAL (2, G1 
DO 70 K=3, LMAXM1 
PVNS=0 . 

DO 75 J=1 , JMAX 
DO 75 1=1, IMAX 

PVNS=PVNS+DUDX(I,J,K)XVa,J,K) 
75 CONTINUE 

PVNS=PVNSXCIJ - 
PVNST (K) =PVNSTCK)+PVNS 
70 CONTINUE 

CALL PARTIAL (3, DIVC) 

DO SO K=3,LMAXM1 
PW=0. 

DO 85 J=1 , JMAX 
DO 85 1=1, IMAX 
PW=PW+DUDX(I,J,K1XU(I,J,K) 

85 CONTINUE 
PW=PWXCIJ 
PWT(K)=PWT(K)+PW 
80 CONTINUE 

CALL PARTIAL(3,G) 

DO 90 K=3,LMAXM1 
PWNS=0. 

DO 95 J=l, JMAX 
DO 95 1=1, IMAX 

PWNS=PWNS+DUDXCI,J,IOXWCI, J,K) 
95 CONTINUE 

pwns=pwns*cij 

PWNST(K) =PWNST(K)+PWNS 
90 CONTINUE 

CALL PARTIAL (1 , U) 

DO 100 K=3, LMAXM1 
PDU=0 . 

PDUN=0 . 

DO 105 J=1 , JMAX 
DO 105 1=1, IMAX 

PDU=PDU+DUDX(I, J,K)XDIVC(I,J,K) 
PDUN=PDUN+DUDX(I,J,K)XG(I,J,K) 
105 CONTINUE 
PDU=PDUXCIJ 
PDUN=PDUNXCIJ 

PDUN=PDUN*CIJ 
PDUT(K)=PDUT(K)+PDU 
PDUNT(K) =PDUNT(K)+PDUN 
100 CONTINUE 

CALL PARTIALC2, V) 

DO 110 K=3, LMAXM1 
PDV=0 . 

PDVN=0. 

DO 115 J=l> JMAX 
DO 115 1=1, IMAX 
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PDV=PDV+DUDXCI,J,K)*DIVCCI,J,10 
PDVN=PDVN+DUDX(I, J>K)XGCI, J,K) 
115 CONTINUE 

PDV=PDV*CIJ 
PDVN=PDVN*CIJ 
PDVTCK)=PDVT(K)+PDV 
PDVNT (K)=PDVNT(K)+PDVN 
110 CONTINUE 

CALL PARTIAL(3,W) 

DO 120 K=3,LMAXM1 
PDWN=0. 

PDW=0. 


DO 125 I=1,IMAX 
DO 125 J=1,JMAX 

PDW=PDW+DUDXCI, J,K)XDIVCCI, J,K) 

PDWN=PDWN+DUDX(I,J,K)*G(I,J,K) 

125 CONTINUE 

PDW=PDW*CIJ 
PDWN=PDWN*CIJ 
PDWT(K)=PDWTCK)+PDW 
PDWNT (K) =PDWNT(K)+PDL*1N 
120 CONTINUE 
RETURN 
END 

XDECK LTPR 

SUBROUTINE LTPR 

CXXXXX THIS SUBROUTINE PRINTS LONG TIME AVERAGES AT DESIGNATED INTERVALS 
CGMMON/RANGE/LMAXM1 , LMAXM2 , LMAXM3 > LMAXM4 , LMAXM5 
C0MM0N/LTA1/USUMC65),UTSUMC65),ST$UM(65),U2SMT(65),V2SMT(65) 

1 ,W2SMT(65) > PVT (65) , PUT(65) ,PUNST (65) , PVNST(65) >PWNST ( 65) , PUT (65) 

2 » TCONT 

COMMON/LTA2/PDUT(65) , PDVTC65) , PDWTC65) ,PDUNT<65) > PDVNT(65) > PDWNT 
1(65) 

C0MM0N/SGTT/SGST(65) , ETEDC 65) , U2STTC65) » V2STT(65) ,W2STT(65) 
l»TSH6Sj TSCNT 

PRINT 10, TCONT, TSHGS,TSCNT 
10 FORMAT (// , 10X,X COUNTERS X,1P3E14.5) 

Fl=l. /TCONT 
F2 = l ,/TSHGS 


F3-1./TSCNT 

DO 20 K=3,LMAXM1 

A1-UTSUM(K)XF1 

A2=U2SMTCK)XF1 

A3"V2SMT ( K)XF1 

A4-W2SMT (K)XF1 

A5~STSUM(K)*F1 

A6=PUT(K)*F1 

A7~PVT(K)XF1 

A8=PWT(K)*F1 

PRINT 30,A1,A2,A3,A4,A5,A6,A7,A8,K 
20 CONTINUE 
PRINT 40 


40 FORMAT (/////) 


DO 50 K=3, LMAXM1 

A1=PUNST(K)XF1 

A2=PVNST(K)XF1 

A3=PWNST(K)XF1 

A4=SGST(K)XF2 

A5=ETEDCK)XF3 

A6=U2STT(K)*F3 

A7=V2STT(K)*F3 

A8=W2STT(K)*F3 

PRINT 30,A1,A2,A3,A4,A5,A6,A7,A8,K 
50 CONTINUE 

30 FORMAT (3X, 1P8E14.5,I5) 

PRINT 40 

DO 60 K=3» LMAXM1 
A1=PDUT(K)XF1 
A2=PDVT(K)XF1 
A3=PDWT(K)XF1 
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A4=PDUNTCiO*Fl 

A5=PDVNTCK)*F1 

A6=PDUNT(IOKFl 

PRINT 70,A1,A2,A3,A4,A5,A6,K 
60 CONTINUE 
70 F0RMAT(1P6E14.5» 15) 

RETURN 
END ' ' 


PR1GINAX PAGE IS 
Qg POOR QUALITY 
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